2024 - deep_analysis.ipynb

{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "# Deep analysis of SPAN peak calling\n",
    "\n",
    "Consider launch SPAN with `--deep-analysis` command line argument."
   ],
   "id": "697668dfbd591db1"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'\n",
    "\n",
    "from io import StringIO\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",
    "import math\n",
    "from itertools import product\n",
    "import glob"
   ],
   "id": "59f961c43bda7f7f",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Plots",
   "id": "cd438387b136ad6d"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def plot_sensitivity(t, name, x='CandidatesN', y='CandidatesAL', logx=True, logy=True):\n",
    "    plt.figure(figsize=(5, 4))\n",
    "    ax = plt.axes()\n",
    "    sns.lineplot(data=t[t['Gap'] > 0], x=x, y=y, hue='Gap', estimator=None,\n",
    "                 palette='tab20',\n",
    "                 hue_order=list(sorted(t['Gap'].unique())),\n",
    "                 sort=False,\n",
    "                 alpha=0.5,\n",
    "                 ax=ax)\n",
    "    sns.lineplot(data=t[t['Gap'] == 0], x=x, y=y, hue='Gap', estimator=None,\n",
    "                 palette='tab20',\n",
    "                 hue_order=list(sorted(t['Gap'].unique())),\n",
    "                 sort=False,\n",
    "                 legend=False,\n",
    "                 ax=ax)\n",
    "\n",
    "    sp = detect_sensitivity_triangle(t[t['Gap'] == 0])\n",
    "    sns.scatterplot(data=t[t.index == sp[0]],\n",
    "                    x=x, y=y, color='green', s=20,\n",
    "                    legend=False, ax=ax)\n",
    "    sns.scatterplot(data=t[t.index == sp[1]],\n",
    "                    x=x, y=y, color='red', s=40,\n",
    "                    legend=False, ax=ax)\n",
    "    sns.scatterplot(data=t[t.index == sp[2]],\n",
    "                    x=x, y=y, color='blue', s=20,\n",
    "                    legend=False, ax=ax)\n",
    "\n",
    "    if logx:\n",
    "        ax.set(xscale='log')\n",
    "    if logy:\n",
    "        ax.set(yscale='log')\n",
    "    ax.set_title(name)\n",
    "    if len(t['Gap'].unique()) > 1:\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    return sp"
   ],
   "id": "500ac0bd8f7265cb",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def plot_projections(t, name, sp, show_text=False):\n",
    "    for y in ['CandidatesN', 'CandidatesAL', 'CandidatesML', 'SignalNoiseRatio', 'SignalControlRatio']:\n",
    "        if y not in t.columns:\n",
    "            continue\n",
    "        print(y)\n",
    "        plt.figure(figsize=(5, 2))\n",
    "        ax = plt.axes()\n",
    "        tgap = t[t['Gap'] > 0]\n",
    "        if len(tgap) > 0:\n",
    "            sns.lineplot(data=tgap, x='SensitivityN', y=y, estimator=None,\n",
    "                         # marker='o', markersize=3,\n",
    "                         hue='Gap',\n",
    "                         hue_order=list(sorted(t['Gap'].unique())),\n",
    "                         palette='tab20',\n",
    "                         alpha=0.5,\n",
    "                         sort=False,\n",
    "                         ax=ax)\n",
    "        tnogap = t[t['Gap'] == 0]\n",
    "        sns.lineplot(data=tnogap, x='SensitivityN', y=y, estimator=None,\n",
    "                     # marker='o', markersize=3,\n",
    "                     hue='Gap',\n",
    "                     hue_order=list(sorted(t['Gap'].unique())),\n",
    "                     palette='tab20',\n",
    "                     alpha=1,\n",
    "                     sort=False,\n",
    "                     legend=False,\n",
    "                     ax=ax)\n",
    "\n",
    "        if sp is not None:\n",
    "            sns.scatterplot(x=[tnogap['SensitivityN'].values[sp[0]]], y=[tnogap[y].values[sp[0]]],\n",
    "                            color='green', s=20,\n",
    "                            legend=False, ax=ax)\n",
    "            sns.scatterplot(x=[tnogap['SensitivityN'].values[sp[1]]], y=[tnogap[y].values[sp[1]]],\n",
    "                            color='red', s=40,\n",
    "                            legend=False, ax=ax)\n",
    "            sns.scatterplot(x=[tnogap['SensitivityN'].values[sp[2]]], y=[tnogap[y].values[sp[2]]],\n",
    "                            color='blue', s=20,\n",
    "                            legend=False, ax=ax)\n",
    "            if show_text and y in ['CandidatesN', 'CandidatesAL', 'CandidatesML']:\n",
    "                ax.text(tnogap['SensitivityN'].values[sp[0]], 10,\n",
    "                        f\"t1={tnogap['Sensitivity'].values[sp[0]]:.2e}\", fontsize=5)\n",
    "                ax.text(tnogap['SensitivityN'].values[sp[1]], 100,\n",
    "                        f\"t1={tnogap['Sensitivity'].values[sp[0]]:.2e}\", fontsize=5)\n",
    "                ax.text(tnogap['SensitivityN'].values[sp[2]], 10,\n",
    "                        f\"t1={tnogap['Sensitivity'].values[sp[0]]:.2e}\", fontsize=5)\n",
    "\n",
    "        ax.set(yscale='log')\n",
    "        ax.set_title(name)\n",
    "        sens = np.array(tnogap['Sensitivity'])\n",
    "        # First and last labels are out of the plot\n",
    "        sens_ticks = np.array([sens[int(t)] if 0 <= t < len(sens) else np.nan for t in ax.get_xticks()[1:-1]])\n",
    "        sens_ticks[-1] = sens.max()\n",
    "        labels = [''] + [f'{t:.1e}' for t in sens_ticks] + ['']\n",
    "        ax.set_xticklabels(labels)\n",
    "        ax.set_xlabel('Sensitivity')\n",
    "        if len(t['Gap'].unique()) > 1:\n",
    "            sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n"
   ],
   "id": "c75dee826dcbdb2",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Area based triangles",
   "id": "cb09bf259be599c2"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def triangle_signed_square(x1, y1, x2, y2, x3, y3):\n",
    "    return 0.5 * (x1 * y2 - x2 * y1 + x2 * y3 - x3 * y2 + x3 * y1 - x1 * y3)\n",
    "\n",
    "\n",
    "def detect_sensitivity_triangle(t):\n",
    "    im1 = int(len(t) * 0.1)\n",
    "    im2 = int(len(t) * 0.9)\n",
    "    ns = np.log1p(t['CandidatesN'].values)\n",
    "    als = np.log1p(t['CandidatesAL'].values)\n",
    "    # print('Limit', 'Im1', im1, 'Im2', im2)\n",
    "    max_area = 0\n",
    "    i1, i2, i3 = -1, -1, -1\n",
    "    for i in range(im1, im2):\n",
    "        i1mab = find_sensitivity_max_area_between(ns, als, 0, i, -1)\n",
    "        i3mab = find_sensitivity_max_area_between(ns, als, i, len(t) - 1, -1)\n",
    "        # print('Checking', 'I', i,  'I1mab', i1mab, 'I3mab', i3mab)\n",
    "        if i1mab[0] == -1 or i3mab[0] == -1:\n",
    "            continue\n",
    "        # we want area to be balanced between parts, geometric mean is better here\n",
    "        area = math.sqrt(i1mab[1] * i3mab[1])\n",
    "        if area > max_area:\n",
    "            max_area = area\n",
    "            i1, i3, i2 = i1mab[0], i3mab[0], i\n",
    "    # print('Found', 'I1', i1, 'I2', i2, 'I3', i3, max_area)\n",
    "    if i3 == -1 or i2 == -1 or i1 == -1:\n",
    "        print(\"Failed to estimate sensitivity triangle\")\n",
    "        return im1, (im1 + im2) / 2, im2\n",
    "\n",
    "    # Update i3, i1 points to be closer to i2 for more accurate pivot_gap estimation\n",
    "    i3mab = find_sensitivity_max_area_between(ns, als, i3, i2, -1)\n",
    "    if i3mab[0] != -1:\n",
    "        i3 = i3mab[0]\n",
    "    i1mab = find_sensitivity_max_area_between(ns, als, i2, i1, -1)\n",
    "    if i1mab[0] != -1:\n",
    "        i1 = i1mab[0]\n",
    "    # print('Result', 'I1', i1, 'I2', i2, 'I3', i3, max_area)\n",
    "    return i1, i2, i3\n",
    "\n",
    "\n",
    "def find_sensitivity_max_area_between(ns, als, start, end, sign):\n",
    "    if start > end:\n",
    "        return find_sensitivity_max_area_between(ns, als, end, start, sign)\n",
    "    max_i = -1\n",
    "    max_area = 0\n",
    "    start_n, start_al = ns[start], als[start]\n",
    "    end_n, end_al = ns[end], als[end]\n",
    "    for i in range(start + 1, end):\n",
    "        n, al = ns[i], als[i]\n",
    "        area = triangle_signed_square(start_n, start_al, n, al, end_n, end_al)\n",
    "        if area * sign > 0:\n",
    "            continue\n",
    "        area = math.fabs(area)\n",
    "        if area > max_area:\n",
    "            max_i, max_area = i, area\n",
    "    return max_i, max_area"
   ],
   "id": "58a2037e88d57ff3",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Bad quality example",
   "id": "360a773b92c71f73"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "t = StringIO('''Sensitivity\tGap\tCandidatesN\tCandidatesAL\n",
    "-375.1737060546875\t0\t1\t1.0\n",
    "-280.0399200217636\t0\t9\t1.6666666666666667\n",
    "-209.029458995094\t0\t92\t1.4891304347826086\n",
    "-156.02530783606127\t0\t453\t1.6291390728476822\n",
    "-116.46155906622265\t0\t1239\t1.883777239709443\n",
    "-86.93009440741795\t0\t2916\t2.0174897119341564\n",
    "-64.8870011209932\t0\t5615\t2.1814781834372217\n",
    "-48.433433130109364\t0\t10588\t2.3860974688326406\n",
    "-36.15203976516999\t0\t16322\t2.576277416983213\n",
    "-26.984871703631814\t0\t24183\t2.750568581234752\n",
    "-20.142246622637973\t0\t34899\t3.1237857818275594\n",
    "-15.03472402844767\t0\t45585\t3.3377426785126687\n",
    "-11.222329407759858\t0\t53435\t3.938748011602882\n",
    "-8.376653744889191\t0\t47325\t5.995456946645536\n",
    "-6.252563564320888\t0\t27576\t13.580395996518712\n",
    "-4.66708453237974\t0\t24350\t17.816632443531827\n",
    "-3.483639599710967\t0\t23656\t20.143599932363884\n",
    "-2.600283919538606\t0\t23517\t21.657481821660927\n",
    "-1.9409230687273291\t0\t23611\t22.653593663970184\n",
    "-1.4487580876885016\t0\t23776\t23.33983849259758\n",
    "-1.0813926788036474\t0\t23999\t23.867577815742322\n",
    "-0.8071810854467272\t0\t24240\t24.26526402640264\n",
    "-0.602502048953911\t0\t24478\t24.615859138818532\n",
    "-0.4497240155135174\t0\t24753\t24.881994101725045\n",
    "-0.33568631091090884\t0\t25075\t25.060019940179462\n",
    "-0.25056544779870277\t0\t25416\t25.208215297450426\n",
    "-0.18702890642218395\t0\t25772\t25.32826323141394\n",
    "-0.1396034933977802\t0\t26168\t25.388833690003057\n",
    "-0.10420386742181362\t0\t26637\t25.404062018996132\n",
    "-0.07778061795863045\t0\t27109\t25.409900770961674\n",
    "-0.05805758154384941\t0\t27699\t25.306509260262104\n",
    "-0.04333576748533296\t0\t28378\t25.138734230742124\n",
    "-0.032347002641238944\t0\t29185\t24.910604762720574\n",
    "-0.024144687877662503\t0\t30126\t24.61840270862378\n",
    "-0.018022255699405163\t0\t31272\t24.22521744691737\n",
    "-0.013452304794349106\t0\t32592\t23.79120643102602\n",
    "-0.01004116839192559\t0\t34071\t23.340700302309884\n",
    "-0.007495002842736571\t0\t36080\t22.680182926829268\n",
    "-0.005594475206471149\t0\t39225\t21.588833652007647\n",
    "-0.004175869374906439\t0\t43757\t20.204241607057156\n",
    "-0.0031169831651252974\t0\t48401\t19.266853990620028\n",
    "-0.002326601523998155\t0\t53579\t18.580413968159167\n",
    "-0.001736639039965088\t0\t63313\t17.234154123165858\n",
    "-0.0012962748988267487\t0\t82491\t15.135032912681384\n",
    "-9.675750542623222E-4\t0\t110976\t13.33552299596309\n",
    "-7.222244961142835E-4\t0\t155013\t11.45052350448027\n",
    "-5.390881260216058E-4\t0\t212609\t10.160722264814753\n",
    "-4.023901282510646E-4\t0\t235773\t10.838098510007507\n",
    "-3.0035500226806847E-4\t0\t244449\t12.541376728888235\n",
    "-2.241931922623222E-4\t0\t292100\t12.481119479630264\n",
    "-1.6734393326970745E-4\t0\t968964\t5.192648024075198\n",
    "-1.2491009079976728E-4\t0\t617855\t10.359380437157586\n",
    "-9.32363096692581E-5\t0\t592402\t12.330579572655056\n",
    "-6.959413274846464E-5\t0\t569118\t14.391375777958174\n",
    "-5.1946964977399616E-5\t0\t549164\t15.138854695500797\n",
    "-3.8774636076238994E-5\t0\t530365\t16.720228521867018\n",
    "-2.8942449351927953E-5\t0\t507602\t18.508398312063388\n",
    "-2.1603436144233477E-5\t0\t479059\t21.157412761267402\n",
    "-1.6125395869679023E-5\t0\t25070\t42.360430793777425\n",
    "-1.2036436713946952E-5\t0\t25070\t42.360430793777425\n",
    "-8.984325714524885E-6\t0\t25070\t42.360470682090146\n",
    "-6.706146550095041E-6\t0\t25070\t42.360470682090146\n",
    "-5.005651284285609E-6\t0\t25070\t42.360470682090146\n",
    "-3.7363550874853966E-6\t0\t25070\t42.360470682090146\n",
    "-2.788917674629893E-6\t0\t25070\t42.360470682090146\n",
    "-2.0817244650849684E-6\t0\t25070\t42.360470682090146\n",
    "-1.553856102657593E-6\t0\t25070\t42.360470682090146\n",
    "-1.159840712957992E-6\t0\t25070\t42.360470682090146\n",
    "-8.657368447014698E-7\t0\t25070\t42.360470682090146\n",
    "-6.46209669914216E-7\t0\t25070\t42.360470682090146\n",
    "-4.823485797634483E-7\t0\t25070\t42.360470682090146\n",
    "-3.60038178368175E-7\t0\t25070\t42.360470682090146\n",
    "-2.68742348005349E-7\t0\t25070\t42.360470682090146\n",
    "-2.0059664210825222E-7\t0\t25070\t42.360470682090146\n",
    "-1.497308225658031E-7\t0\t25070\t42.360470682090146\n",
    "-1.1176318302543429E-7\t0\t25070\t42.360470682090146\n",
    "-8.34230979696061E-8\t0\t25070\t42.360470682090146\n",
    "-6.226928301838606E-8\t0\t25070\t42.360470682090146\n",
    "-4.6479496710090816E-8\t0\t25070\t42.360470682090146\n",
    "-3.469356815599526E-8\t0\t25070\t42.360470682090146\n",
    "-2.589622858660116E-8\t0\t25070\t42.360470682090146\n",
    "-1.9329653611705986E-8\t0\t25070\t42.360470682090146\n",
    "-1.4428182370225825E-8\t0\t25070\t42.360470682090146\n",
    "-1.0769590117353544E-8\t0\t25070\t42.360470682090146\n",
    "-8.03871675029177E-9\t0\t25070\t42.360470682090146\n",
    "-6.000318144633462E-9\t0\t25070\t42.360470682090146\n",
    "-4.478801648970997E-9\t0\t25070\t42.360470682090146\n",
    "-3.3431001035780287E-9\t0\t25070\t42.360470682090146\n",
    "-2.4953813940635528E-9\t0\t25070\t42.360470682090146\n",
    "-1.8626209532804875E-9\t0\t25070\t42.360470682090146\n",
    "-1.3903112461497955E-9\t0\t25070\t42.360470682090146\n",
    "-1.0377663570068893E-9\t0\t25070\t42.360470682090146\n",
    "-7.746172051170474E-10\t0\t25070\t42.360470682090146\n",
    "-5.781954776352038E-10\t0\t25070\t42.360470682090146\n",
    "-4.315809255841242E-10\t0\t25070\t42.360470682090146\n",
    "-3.22143812140927E-10\t0\t25070\t42.360470682090146\n",
    "-2.404569561554013E-10\t0\t25070\t42.360470682090146\n",
    "-1.7948365166246456E-10\t0\t25070\t42.360470682090146\n",
    "-1.339715087854371E-10\t0\t25070\t42.360470682090146\n",
    "-9.999999999999925E-11\t0\t25070\t42.360470682090146''')\n",
    "t = pd.read_csv(t, sep='\\t')\n",
    "t['SensitivityN'] = range(len(t))\n",
    "# For log\n",
    "t['CandidatesN'] += 1\n",
    "t['CandidatesAL'] += 1\n",
    "t.sample(3)"
   ],
   "id": "a767f97f724d84fe",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sp = plot_sensitivity(t, 'Mix noise bad quality')\n",
    "display(sp)"
   ],
   "id": "db85679a0993084c",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "plot_projections(t, 'Mix noise bad quality', sp)",
   "id": "382374f7b233ab3c",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Peak infos",
   "id": "4c05e3126dba421b"
  },
  {
   "metadata": {},
   "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_ABF = ['k4me3', 'k27ac', 'k4me1', 'k27me3', 'k36me3']\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",
    "IMMGEN_PATH = os.path.expanduser('~/data/2024_Immgen')\n",
    "IMMGEN_CELLS = [\n",
    "\n",
    "    'Thymus_Thymic_epithelial_Cell', 'Thymus_Double_Negative_Thymocytes',\n",
    "                'Thymus_Immature_Single_Positive_Thymocytes',\n",
    "                'Thymus_CD4_SP_Thymocytes_ATAC',\n",
    "                'Thymus_CD8_SP_Thymocytes',\n",
    "                'Spleen_CD4_Naive_T_Cell',\n",
    "                'Spleen_CD8_Naive_T_Cell',\n",
    "                'Spleen_Activated_T_Cell',\n",
    "                'Colon_CD4_Treg_Cell',\n",
    "                'Bone_marrow_Neutrophil',\n",
    "                'Spleen_Neutrophil_ATAC',\n",
    "                'Peritoneal_cavity_Neutrophil',\n",
    "                'Spleen_NK_Cell',\n",
    "                'Bone_marrow_NK_Cell',\n",
    "                'Thymus_gdT_Cell',\n",
    "                'Lymph_nodes_gdT_Cell',\n",
    "                'Spleen_CD8_T_Cell',\n",
    "                '']  ## Part, longest first\n",
    "IMMGEN_REPS = ['']\n",
    "\n",
    "CTCF_PATH = os.path.expanduser('~/data/2024_TFs')\n",
    "CTCF_CELLS = ['H1']  ## Longest first\n",
    "CTCF_REPS = ['']\n",
    "\n",
    "Y20O20_PATH = os.path.expanduser('~/data/2018_chipseq_y20o20')\n",
    "Y20O20_CELLS = ['']\n",
    "Y20O20_REPS = [f'OD{i}' for i in range(1, 30)] + [f'YD{i}' for i in range(1, 30)]"
   ],
   "id": "e1507a4d0d12b365",
   "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'"
   ],
   "id": "bb22c251363db6e5",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def load_peaks_infos(path, suffix, modifications, cells, replicates):\n",
    "    ts = []\n",
    "    vals = {}\n",
    "    for f in os.listdir(path):\n",
    "        if not f.endswith(f'{suffix}.txt'):\n",
    "            continue\n",
    "        if 'ATAC_seq' in f:\n",
    "            rep = re.sub('_.*', '', os.path.basename(f))\n",
    "            cell = re.sub('(SRR[0-9]+_)|(_ATAC.*)', '', os.path.basename(f))\n",
    "        else:\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.lower() in f.lower()), None)\n",
    "        if mod and cell is not None and rep is not None:\n",
    "            txt_path = os.path.join(path, f)\n",
    "            print(txt_path, mod, cell, rep)\n",
    "            with open(txt_path, 'r') as file:\n",
    "                txt = file.read().rstrip()\n",
    "                vals = {}\n",
    "                for (c, text) in [\n",
    "                    ('peaks', 'Count'),\n",
    "                    ('length_mean', 'Mean length'),\n",
    "                    ('length_median', 'Median length'),\n",
    "                    ('model_signal', 'Signal mean'),\n",
    "                    ('model_noise', 'Noise mean'),\n",
    "                    ('model_signal_to_noise', 'Signal to noise'),\n",
    "                    ('frip', 'FRIP'),\n",
    "                    ('log_null_pvals_mean', 'LogNullPVals mean'),\n",
    "                    ('log_null_pvals_std', 'LogNullPVals std'),\n",
    "                    ('coverage', 'Treatment coverage'),\n",
    "                    ('control_coverage', 'Control coverage'),\n",
    "                    ('control_scale', 'Control scale'),\n",
    "                    ('beta', 'Beta'),\n",
    "                    ('min_correlation', 'Min control correlation'),\n",
    "                    ('coverage_non_zero', 'Coverage >0 %'),\n",
    "                    ('coverage_max', 'Coverage >0 max'),\n",
    "                    ('coverage_mean', 'Coverage >0 mean'),\n",
    "                    ('coverage_median', 'Coverage >0 median'),\n",
    "                    ('coverage_std', 'Coverage >0 std'),\n",
    "                    ('roughness', 'Track roughness'),\n",
    "                    ('autocorrelation_average_score', 'Average autocorrelation score'),\n",
    "                    ('sensitivity_before_merge', 'Sensitivity beforeMerge'),\n",
    "                    ('sensitivity_before_merge_idx', 'Sensitivity beforeMerge index'),\n",
    "                    ('sensitivity_stable', 'Sensitivity stable'),\n",
    "                    ('sensitivity_stable_idx', 'Sensitivity stable index'),\n",
    "                    ('sensitivity_before_noise', 'Sensitivity beforeNoise'),\n",
    "                    ('sensitivity_before_noise_idx', 'Sensitivity beforeNoise index'),\n",
    "                    ('minimal_additional', 'Minimal additional'),\n",
    "                    ('minimal_additional_idx', 'Minimal additional index'),\n",
    "                    ('fragmentation_average_score', 'Average fragmentation score'),\n",
    "                    ('signal_density', 'Candidates signal density'),\n",
    "                    ('noise_density', 'Candidates noise density'),\n",
    "                    ('signal_to_noise', 'Coverage signal to noise'),\n",
    "                    ('signal_to_control', 'Coverage signal to control'),\n",
    "                ]:\n",
    "                    try:\n",
    "                        if text in txt:\n",
    "                            x = re.sub('\\n(\\n|.)*', '', re.sub(f'(.|\\n)*{text}: ', '', txt)).replace(',', '')\n",
    "                            vals[c] = float(x)\n",
    "                    except Exception as e:\n",
    "                        print(f'Failed to process {c}: {text}: {x}', e)\n",
    "                fit_snr_adjusted = re.sub('\\n(\\n|.)*', '',\n",
    "                                          re.sub('(.|\\n)*Out of low noise level down: ', '', txt)) == 'true'\n",
    "                vals['fit_snr_adjusted'] = fit_snr_adjusted\n",
    "                fit_low_adjusted = re.sub('\\n(\\n|.)*', '',\n",
    "                                          re.sub('(.|\\n)*Out of signal-to-noise range down: ', '', txt)) == 'true'\n",
    "                vals['fit_low_adjusted'] = fit_low_adjusted\n",
    "\n",
    "                ts.append((txt_path, mod, cell, rep, *vals.values()))\n",
    "    df = pd.DataFrame(ts, columns=['file', 'modification', 'cell', 'replicate', *vals.keys()])\n",
    "    return df"
   ],
   "id": "2029ebf68151c1ce",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "gse26320_infos = load_peaks_infos(os.path.join(GSE26320_PATH, 'span'), '.peak',\n",
    "                                  MODIFICATIONS, GSE26320_CELLS, GSE26320_REPS)\n",
    "gse26320_infos['dataset'] = 'ENCODE'\n",
    "immune_infos = load_peaks_infos(os.path.join(IMMUNE_PATH, 'span'), '.peak',\n",
    "                                MODIFICATIONS, IMMUNE_CELLS, IMMUNE_REPS)\n",
    "immune_infos['dataset'] = 'Roadmap'\n",
    "immgen_infos = load_peaks_infos(os.path.join(IMMGEN_PATH, 'span'), '.peak',\n",
    "                                ['ATAC'], IMMGEN_CELLS, IMMGEN_REPS)\n",
    "immgen_infos['dataset'] = 'Immgen'\n",
    "ctcf_infos = load_peaks_infos(os.path.join(CTCF_PATH, 'span'), '.peak',\n",
    "                                ['CTCF'], CTCF_CELLS, CTCF_REPS)\n",
    "ctcf_infos['dataset'] = 'ENCODE'\n",
    "\n",
    "y20o20_infos = load_peaks_infos(os.path.join(Y20O20_PATH, 'span'), '.peak',\n",
    "                                MODIFICATIONS_ABF, Y20O20_CELLS, Y20O20_REPS)\n",
    "y20o20_infos['dataset'] = 'ABF'\n",
    "\n",
    "full_infos = pd.concat([gse26320_infos, immune_infos, immgen_infos, ctcf_infos, y20o20_infos]).reset_index(drop=True)\n",
    "update_abf_modifications(full_infos)\n",
    "# full_infos = pd.concat([gse26320_infos, immune_infos]).reset_index(drop=True)\n",
    "full_infos.sample(3)"
   ],
   "id": "c473a0c23ba75413",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "full_infos['modification'].unique()",
   "id": "71f98535d980a4f",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Remove outliers\n",
    "OUTLIERS = {\n",
    "    'H3K4me3:': {'NK': ['']},\n",
    "    'H3K27ac': {'': ['YD1', 'YD6']},\n",
    "    'H3K27me3': {\n",
    "        'TCell':[''],\n",
    "        'BCell': [''],\n",
    "        '': ['OD9', 'YD1', 'YD2', 'YD3', 'YD4', 'YD9', 'YD10', 'YD11']},\n",
    "    'H3K36me3': {'': ['OD3', 'OD6', 'OD12', 'OD18', 'OD20', 'YD1', 'YD3', 'YD4', 'YD5']},\n",
    "    'H3K4me3': {'': ['OD6', 'OD7', 'OD14', 'YD2', 'YD3', 'YD10', 'YD14']}\n",
    "}\n",
    "full_infos['outlier'] = [\n",
    "    any(m == mo and c == co and r == ro\n",
    "           for mo, crlo in OUTLIERS.items()\n",
    "                for co, rlo in crlo.items()\n",
    "                    for ro in rlo) for m, c, r in zip(full_infos['modification'], full_infos['cell'], full_infos['replicate'])\n",
    "]\n",
    "print('Total outliers', sum(full_infos['outlier']))\n",
    "df_peaks = full_infos[~full_infos['outlier']].copy()\n",
    "print('Good peaks', len(df_peaks))"
   ],
   "id": "41137ed160506b2f",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sns.set_style('whitegrid')\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "sns.barplot(data=df_peaks, x='modification', y='peaks', hue='dataset',\n",
    "            capsize=.1, errwidth=2, edgecolor=\"black\",\n",
    "            # order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "            order=MODIFICATIONS,\n",
    "            # hue_order=['ENCODE', 'Roadmap', 'ABF', 'Immgen'],\n",
    "            hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "            ax=ax)\n",
    "sns.swarmplot(data=df_peaks, x='modification', y='peaks', hue='dataset',\n",
    "              dodge=True, size=2, palette='dark:black', alpha=0.5,\n",
    "              # order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "              order=MODIFICATIONS,\n",
    "              # hue_order=['ENCODE', 'Roadmap', 'ABF', 'Immgen'],\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('Peaks')\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()"
   ],
   "id": "ff2faba4e27a3e53",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for out in ['fit_snr_adjusted', 'fit_low_adjusted']:\n",
    "    print(out)\n",
    "    t = full_infos[['dataset', 'modification']].copy()\n",
    "    t['number'] = [1 if o == True else 0 for o in full_infos[out]]\n",
    "    t1 = t.groupby(['dataset', 'modification']).sum().reset_index()\n",
    "    t2 = t.groupby(['dataset', 'modification']).count().reset_index()\n",
    "    t1['number'] /= t2['number']\n",
    "    t1['number'] *= 100\n",
    "    plt.figure(figsize=(3, 1))\n",
    "    ax = plt.axes()\n",
    "    sns.barplot(data=t1, x='modification', y='number',\n",
    "                hue='dataset',\n",
    "                order=['CTCF', 'ATAC' ] + MODIFICATIONS,\n",
    "                hue_order=['ENCODE', 'Roadmap', 'Immgen'],\n",
    "                ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.set_ylabel('%')\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    # plt.tight_layout()\n",
    "    ax.set_title(f'{out}')\n",
    "    plt.show()"
   ],
   "id": "98afce65fce8aa11",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Autocorrelation vs fragmentation",
   "id": "9971e7dfd88302ff"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(2.5, 3))\n",
    "ax = plt.axes()\n",
    "sns.barplot(data=df_peaks, x='modification', y='autocorrelation_average_score',\n",
    "            hue='modification', legend=False,\n",
    "            order=['CTCF', 'ATAC' ] + MODIFICATIONS,\n",
    "            capsize=.2, err_kws={'linewidth': 2},\n",
    "            ax=ax)\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "ax.set_title('Average PEP autocorrelation')\n",
    "ax.set_ylabel('Autocorrelation')\n",
    "ax.set_ylabel('Experiment')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "d05c4c141b7a3b4d",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for y in ['autocorrelation_average_score', 'fragmentation_average_score']:\n",
    "    plt.figure(figsize=(5, 2.5))\n",
    "    ax = plt.axes()\n",
    "    sns.boxplot(df_peaks,\n",
    "                x='modification', y=y,\n",
    "                hue='dataset',\n",
    "                order=['CTCF', 'ATAC' ] + MODIFICATIONS,\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF', 'Immgen'],\n",
    "                ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.set_title(y)\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "b00a704e938fd0d2",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 3.5))\n",
    "ax = plt.axes()\n",
    "sns.scatterplot(df_peaks, x='autocorrelation_average_score', y='length_mean',\n",
    "                hue='modification',\n",
    "                hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "                style='dataset',\n",
    "                alpha=0.8,\n",
    "                ax=ax)\n",
    "ax.set_xlabel('autocorrelation')\n",
    "ax.set_ylabel('average length')\n",
    "ax.set_title('Autocorrelation vs length')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "31f4251a232b0e95",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 3.5))\n",
    "ax = plt.axes()\n",
    "sns.kdeplot(\n",
    "    df_peaks, x='autocorrelation_average_score', y='fragmentation_average_score',\n",
    "    hue='modification',\n",
    "    hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "    alpha=0.2,\n",
    "    common_norm=False, common_grid=False,\n",
    "    thresh=.1,\n",
    "    fill=True,\n",
    "    warn_singular=False,\n",
    "    ax = ax,\n",
    ")\n",
    "sns.scatterplot(df_peaks, x='autocorrelation_average_score', y='fragmentation_average_score',\n",
    "                hue='modification',\n",
    "                hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "                style='dataset',\n",
    "                alpha=0.8,\n",
    "                ax=ax)\n",
    "ax.set_xlabel('narrow ← AC ~ width → broad')\n",
    "ax.set_ylabel('high ← fragmentation → low')\n",
    "ax.set_title('Autocorrelation vs Fragmentation')\n",
    "ax.set_xlim(-0.1, 1.1)\n",
    "ax.set_ylim(-0.1, 1.1)\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "55bfd8dafbfa6d98",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "\n",
    "* Empirical almost no fragmentation threshold<br>\n",
    "    `SPAN_DEFAULT_FRAGMENTATION_MAX_THRESHOLD = 0.8`\n",
    "\n",
    "* Additional gap for tracks with high fragmentation<br>\n",
    "    `SPAN_DEFAULT_FRAGMENTATION_COMPENSATION_GAP = 50`  // x default bin 100bp\n",
    "\n",
    "```\n",
    "        if (avgFragmentation < SPAN_DEFAULT_FRAGMENTATION_MAX_THRESHOLD) {\n",
    "            val fragmentationGap =\n",
    "                ceil((SPAN_DEFAULT_FRAGMENTATION_MAX_THRESHOLD - avgFragmentation) /\n",
    "                        SPAN_DEFAULT_FRAGMENTATION_MAX_THRESHOLD * SPAN_DEFAULT_FRAGMENTATION_COMPENSATION_GAP).toInt()\n",
    "        }\n",
    "```"
   ],
   "id": "26b0e5112a07ca33"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 3.5))\n",
    "ax = plt.axes()\n",
    "sns.scatterplot(df_peaks, x='autocorrelation_average_score', y='fragmentation_average_score',\n",
    "                hue='peaks',\n",
    "                style='dataset',\n",
    "                alpha=0.8,\n",
    "                ax=ax)\n",
    "ax.set_xlabel('narrow ← AC ~ width → broad')\n",
    "ax.set_ylabel('high ← fragmentation → low')\n",
    "ax.set_title('Peaks number')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "920c35ce3ec9ba30",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 3.5))\n",
    "ax = plt.axes()\n",
    "sns.scatterplot(df_peaks, x='autocorrelation_average_score', y='fragmentation_average_score',\n",
    "                hue='length_mean',\n",
    "                style='dataset',\n",
    "                alpha=0.8,\n",
    "                ax=ax)\n",
    "ax.set_xlabel('narrow ← AC ~ width → broad')\n",
    "ax.set_ylabel('high ← fragmentation → low')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.set_title('Average length')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "188e9a697842755",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "47aeae39acdc06f8",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# PEP threshold sensitivity",
   "id": "7b801d30ec9be1ea"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def load_peaks_sensitivity(path, modifications, cells, replicates):\n",
    "    ts = []\n",
    "    # print(os.listdir(path))\n",
    "    for f in tqdm(os.listdir(path)):\n",
    "        if not f.endswith('.peak.sensitivity.tsv'):\n",
    "            continue\n",
    "        # print(f)\n",
    "        if 'ATAC_seq' in f:\n",
    "            cell = re.sub('(SRR[0-9]+_)|(_ATAC.*)', '', os.path.basename(f))\n",
    "            rep = re.sub('_.*', '', os.path.basename(f))\n",
    "        else:\n",
    "            cell = 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), None)\n",
    "        mod = next((m for m in modifications if m.lower() in f.lower()), None)\n",
    "        # print(f, cell, mod, rep)\n",
    "        if cell and mod and rep is not None:\n",
    "            try:\n",
    "                tsv_path = os.path.join(path, f)\n",
    "                print(tsv_path, mod, cell, rep)\n",
    "                t = pd.read_csv(tsv_path, sep='\\t')\n",
    "                t['file'] = tsv_path\n",
    "                t['modification'] = mod\n",
    "                t['cell'] = cell\n",
    "                t['replicate'] = rep\n",
    "                for g, tg in t.groupby('Gap'):\n",
    "                    tg = tg.sort_values(by=['Sensitivity']).copy().reset_index(drop=True)\n",
    "                    tg.drop_duplicates(inplace=True)\n",
    "                    tg['SensitivityN'] = range(len(tg))\n",
    "                    ts.append(tg)\n",
    "            except Exception as e:\n",
    "                print(e)\n",
    "    df = pd.concat(ts).reset_index(drop=True)\n",
    "    return df\n"
   ],
   "id": "5604d9abe0710442",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "gse26320_sensitivity = load_peaks_sensitivity(GSE26320_PATH + '/span', MODIFICATIONS, GSE26320_CELLS, GSE26320_REPS)\n",
    "gse26320_sensitivity['dataset'] = 'ENCODE'\n",
    "print('ENCODE')\n",
    "# display(gse26320_sensitivity.head())\n",
    "\n",
    "immune_sensitivity = load_peaks_sensitivity(IMMUNE_PATH + '/span', MODIFICATIONS, IMMUNE_CELLS, IMMUNE_REPS)\n",
    "immune_sensitivity['dataset'] = 'Roadmap'\n",
    "print('Immune')\n",
    "# display(immune_sensitivity.head())\n",
    "\n",
    "immgen_sensitivity = load_peaks_sensitivity(IMMGEN_PATH + '/span', ['ATAC'], IMMGEN_CELLS, IMMGEN_REPS)\n",
    "immgen_sensitivity['dataset'] = 'ImmGen'\n",
    "print('ImmGen')\n",
    "# display(immgen_sensitivity.head())\n",
    "\n",
    "ctcf_sensitivity = load_peaks_sensitivity(CTCF_PATH + '/span', ['CTCF'], CTCF_CELLS, CTCF_REPS)\n",
    "ctcf_sensitivity['dataset'] = 'CTCF'\n",
    "print('CTCF')\n",
    "# display(ctcf_sensitivity.head())\n",
    "\n",
    "full_sensitivity = pd.concat([gse26320_sensitivity, immune_sensitivity, immgen_sensitivity, ctcf_sensitivity]).reset_index(drop=True)\n",
    "# For log scale\n",
    "full_sensitivity.fillna(0, inplace=True)\n",
    "full_sensitivity['CandidatesN'] += 1\n",
    "full_sensitivity['CandidatesAL'] += 1\n",
    "full_sensitivity.sample(3)"
   ],
   "id": "919489a2eebc1af8",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Examples",
   "id": "ece4d7b42f487d5c"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r, ds) in EXAMPLES:\n",
    "    t = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                         (full_sensitivity['modification'] == m) &\n",
    "                         (full_sensitivity['cell'] == c) &\n",
    "                         (full_sensitivity['replicate'] == r)].copy().reset_index(drop=True)\n",
    "    if len(t) == 0:\n",
    "        continue\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    # plot_sensitivity(t, name, 'CandidatesAL', 'CandidatesML', logx=True, logy=True)\n",
    "    # plot_sensitivity(t, name, 'CandidatesN', 'SignalNoiseRatio', logx=True, logy=False)\n",
    "    # plot_sensitivity(t, name, 'CandidatesN', 'SignalControlRatio', logx=True, logy=False)\n",
    "    # plot_sensitivity(t, name, 'SignalControlRatio', 'SignalNoiseRatio', logx=False, logy=False)\n",
    "    sp = plot_sensitivity(t, name)\n",
    "    plot_projections(t, name, sp)"
   ],
   "id": "724c14668fa2b453",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r, ds) in [\n",
    "    ('H3K27me3', 'Monocyte', '', 'Roadmap'),\n",
    "]:\n",
    "    t = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                         (full_sensitivity['modification'] == m) &\n",
    "                         (full_sensitivity['cell'] == c) &\n",
    "                         (full_sensitivity['replicate'] == r) &\n",
    "                         (full_sensitivity['Gap'] == 0)].copy().reset_index(drop=True)\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    plt.figure(figsize=(5, 2))\n",
    "    ax = plt.axes()\n",
    "    sns.lineplot(data=t, x='SensitivityN', y='CandidatesN',\n",
    "             sort=False,\n",
    "             estimator=None,\n",
    "             ax=ax)\n",
    "    ax.set(yscale='log')\n",
    "    sens = np.array(t['Sensitivity'])\n",
    "    # First and last labels are out of the plot\n",
    "    sens_ticks = np.array([sens[int(t)] if 0 <= t < len(sens) else np.nan for t in ax.get_xticks()[1:-1]])\n",
    "    sens_ticks[-1] = sens.max()\n",
    "    labels = [''] + [f'{t:.1e}' for t in sens_ticks] + ['']\n",
    "    ax.set_xticklabels(labels)\n",
    "    ax.set_title('Number of candidates')\n",
    "    ax.set_xlabel('PEP threshold')\n",
    "    ax.set_ylabel('Number')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n"
   ],
   "id": "403c7753e204f2ef",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Per modification",
   "id": "25ea34520e5b5b3b"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "dfg = full_sensitivity[(full_sensitivity['Gap'] == 0)].copy()\n",
    "dfg.sort_values(by=['dataset', 'modification', 'cell', 'replicate', 'SensitivityN'], inplace=True)\n",
    "plt.figure(figsize=(5, 4))\n",
    "ax = plt.axes()\n",
    "for i, (_, t) in enumerate(dfg.groupby(['dataset', 'modification', 'cell', 'replicate'])):\n",
    "    sns.lineplot(data=t, x='CandidatesN', y='CandidatesAL',\n",
    "                 hue='modification',\n",
    "                 hue_order=['ATAC'] + MODIFICATIONS,\n",
    "                 sort=False,\n",
    "                 estimator=None,\n",
    "                 # errorbar=('se', .95),\n",
    "                 alpha=0.3,\n",
    "                 legend = i == 0,\n",
    "                 ax=ax)\n",
    "ax.set(xscale='log')\n",
    "ax.set(yscale='log')\n",
    "ax.set_title('Sensitivity')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "55e926de535d6683",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# dfg = full_sensitivity[full_sensitivity['Gap'] == 0]\n",
    "# for y in ['CandidatesN', 'CandidatesAL', 'SignalNoiseRatio', 'SignalControlRatio']:\n",
    "#     print(y)\n",
    "#     plt.figure(figsize=(8, 2))\n",
    "#     ax = plt.axes()\n",
    "#     sns.lineplot(data=dfg, x='Sensitivity', y=y, hue='modification',\n",
    "#                  hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "#                  alpha=0.5,\n",
    "#                  sort=False,\n",
    "#                  errorbar=('se', .95),\n",
    "#                  ax=ax)\n",
    "#     ax.set(xscale='log')\n",
    "#     ax.set(yscale='log')\n",
    "#     ax.invert_xaxis()\n",
    "#     ax.set_title(y)\n",
    "#     sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "#     plt.tight_layout()\n",
    "#     plt.show()"
   ],
   "id": "f95be58131e59826",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# plt.rc('font', size=5)\n",
    "# dfg = full_sensitivity[full_sensitivity['Gap'] == 0]\n",
    "# for y in ['CandidatesN', 'CandidatesAL', 'SignalNoiseRatio', 'SignalControlRatio']:\n",
    "#     print(y)\n",
    "#     plt.figure(figsize=(3, 6))\n",
    "#     axs = [plt.subplot(7, 1, i + 1) for i in range(7)]\n",
    "#     for i, m in enumerate(['CTCF', 'ATAC'] + MODIFICATIONS):\n",
    "#         ax = axs[i]\n",
    "#         sns.lineplot(data=dfg[dfg['modification'] == m], x='Sensitivity', y=y, hue='dataset',\n",
    "#                      alpha=0.5,\n",
    "#                      sort=False,\n",
    "#                      errorbar=('se', .95),\n",
    "#                      ax=ax)\n",
    "#         ax.set(xscale='log')\n",
    "#         ax.set(yscale='log')\n",
    "#         ax.invert_xaxis()\n",
    "#         ax.set_title(f'{m} {y}')\n",
    "#         sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "#     plt.tight_layout()\n",
    "#     plt.show()\n",
    "# plt.rc('font', size=8)"
   ],
   "id": "db092ba3ed0d2cbe",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Individual",
   "id": "c526d7cec0080765"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# import plotly.express as px\n",
    "#\n",
    "# full_sensitivity['experiment'] = full_sensitivity['dataset'] + ' ' + full_sensitivity['cell'] + ' ' + full_sensitivity['replicate']\n",
    "# for m in ['ATAC'] + MODIFICATIONS:\n",
    "#     dfm = full_sensitivity[(full_sensitivity['modification'] == m) &\n",
    "#                            (full_sensitivity['Gap'] == 0)]\n",
    "#     if len(dfm) == 0:\n",
    "#         continue\n",
    "#     fig = px.line(dfm, x='CandidatesN', y='CandidatesAL', color='experiment',\n",
    "#                   hover_data=['experiment', 'Sensitivity', 'CandidatesN', 'CandidatesAL'],\n",
    "#                   log_x=True, log_y=True,\n",
    "#                   title=f'{m}',\n",
    "#                   width=800, height=600)\n",
    "#     fig.show()"
   ],
   "id": "99110d59ce266517",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# import plotly.express as px\n",
    "#\n",
    "# full_sensitivity['experiment'] = full_sensitivity['dataset'] + ' ' + full_sensitivity['cell'] + ' ' + full_sensitivity['replicate']\n",
    "# plt.rc('font', size=4)\n",
    "# for y in ['CandidatesN', 'CandidatesAL', 'SignalNoiseRatio']:\n",
    "#     print(y)\n",
    "#     for m in ['ATAC'] + MODIFICATIONS:\n",
    "#         dfm = full_sensitivity[(full_sensitivity['modification'] == m) &\n",
    "#                                (full_sensitivity['Gap'] == 0)]\n",
    "#         if len(dfm) == 0:\n",
    "#             continue\n",
    "#         fig = px.line(dfm, x='Sensitivity', y=y, color='experiment',\n",
    "#                       hover_data=['experiment', 'Sensitivity', y],\n",
    "#                       log_x=True, log_y=True,\n",
    "#                       title=f'{m} {y}',\n",
    "#                       width=1000, height=400)\n",
    "#         fig['layout']['xaxis']['autorange'] = \"reversed\"\n",
    "#         fig.show()\n",
    "# plt.rc('font', size=8)"
   ],
   "id": "65f74cded62de399",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "5d799285f3b1e031",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Triangles",
   "id": "e48f82a862338003"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sps = {}\n",
    "for m in ['CTCF', 'ATAC'] + MODIFICATIONS:\n",
    "    for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "        dfdsm = full_sensitivity[(full_sensitivity['dataset'] == ds) & (full_sensitivity['modification'] == m)]\n",
    "        for i, ((cell, rep), t) in enumerate(dfdsm.groupby(['cell', 'replicate'])):\n",
    "            print(ds, m, cell, rep)\n",
    "            tt = t[t['Gap'] == 0].sort_values(by=['Sensitivity']).copy().reset_index(drop=True)\n",
    "            sp = detect_sensitivity_triangle(tt)\n",
    "            name = f'{ds} {m} {cell} {rep}'\n",
    "            sps[name] = sp"
   ],
   "id": "88d2fad0f757f5c4",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## All",
   "id": "a4205006ba0e012d"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for m in ['CTCF', 'ATAC'] + MODIFICATIONS:\n",
    "    for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "        print(ds, m)\n",
    "        dfdsm = full_sensitivity[(full_sensitivity['dataset'] == ds) & (full_sensitivity['modification'] == m)]\n",
    "        t = dfdsm\n",
    "        minx, maxx = t['CandidatesN'].min(), t['CandidatesN'].max()\n",
    "        minx, maxx = 1, maxx + (maxx - minx) * 0.1\n",
    "        miny, maxy = t['CandidatesAL'].min(), t['CandidatesAL'].max()\n",
    "        miny, maxy = 1, maxy + (maxy - miny) * 0.1\n",
    "        n = min(len(dfdsm[['cell', 'replicate']].copy().drop_duplicates()), 18)\n",
    "        print(n)\n",
    "        if n <= 6:\n",
    "            ncols, nrows = n, 1\n",
    "        else:\n",
    "            ncols = int(math.floor(math.sqrt(n) + 2))\n",
    "            nrows = int(math.ceil(n / ncols))\n",
    "        plt.figure(figsize=(3 * ncols + 1, 3 * nrows))\n",
    "        axs = [plt.subplot(nrows, ncols, i + 1) for i in range(n)]\n",
    "        for i, ((cell, rep), t) in enumerate(dfdsm.groupby(['cell', 'replicate'])):\n",
    "            if i == 18:\n",
    "                break\n",
    "            tt = t\n",
    "            ax = axs[i]\n",
    "            sns.lineplot(data=tt, x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap',\n",
    "                         palette='tab20',\n",
    "                         sort=False,\n",
    "                         alpha=0.5,\n",
    "                         ax=ax)\n",
    "\n",
    "            tt = t[t['Gap'] == 0].sort_values(by=['Sensitivity']).copy().reset_index(drop=True)\n",
    "            sp = detect_sensitivity_triangle(tt)\n",
    "            st = full_infos[(full_infos['dataset'] == ds) &\n",
    "                            (full_infos['modification'] == m) &\n",
    "                            (full_infos['cell'] == cell) &\n",
    "                            (full_infos['replicate'] == rep)]\n",
    "            name = f'{ds} {m} {cell} {rep}'\n",
    "            sp = sps[name] if name in sps else None\n",
    "            if sp is not None:\n",
    "                sns.scatterplot(data=tt[tt.index.isin([sp[0],])],\n",
    "                                x=\"CandidatesN\", y=\"CandidatesAL\", color='green', s=20,\n",
    "                                legend=False, ax=ax)\n",
    "                sns.scatterplot(data=tt[tt.index.isin([sp[1],])],\n",
    "                                x=\"CandidatesN\", y=\"CandidatesAL\", color='red', s=40,\n",
    "                                legend=False, ax=ax)\n",
    "                sns.scatterplot(data=tt[tt.index.isin([sp[2],])],\n",
    "                                x=\"CandidatesN\", y=\"CandidatesAL\", color='blue', s=20,\n",
    "                                legend=False, ax=ax)\n",
    "            ax.set(yscale='log')\n",
    "            ax.set(xscale='log')\n",
    "            ax.set_title(f'{m} {cell} {rep}')\n",
    "            ax.set_xlim(minx, maxx)\n",
    "            ax.set_ylim(miny, maxy)\n",
    "            if i == n - 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()\n"
   ],
   "id": "94590c32c3f1dd12",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Figures",
   "id": "8ba4eab689ef09a6"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "EXAMPLES = [\n",
    "    ('CTCF', 'H1', ''),\n",
    "    ('ATAC', 'Blood_Monocyte', 'SRR5799491', 'ImmGen'),\n",
    "    ('H3K4me1', 'CD34', '', 'Roadmap'),\n",
    "    ('H3K4me3', 'K562', 'rep1', 'ENCODE'),\n",
    "    ('H3K27ac', 'K562', 'rep1', 'ENCODE'),\n",
    "    ('H3K27me3', 'Monocyte', '', 'Roadmap'),\n",
    "    ('H3K36me3', 'CD4', '', 'Roadmap'),\n",
    "]"
   ],
   "id": "6e0e237fd78f1b73",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for (m, c, r, ds) in EXAMPLES:\n",
    "    t = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                         (full_sensitivity['modification'] == m) &\n",
    "                         (full_sensitivity['cell'] == c) &\n",
    "                         (full_sensitivity['replicate'] == r) &\n",
    "                         (full_sensitivity['Gap'] == 0)].copy().reset_index(drop=True)\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    sp = sps[name] if name in sps else None\n",
    "    if sp is not None:\n",
    "        t['S2CN0'] = t[t.index == sp[0]]['CandidatesN'].values[0]\n",
    "        t['S2CAL0'] = t[t.index == sp[0]]['CandidatesAL'].values[0]\n",
    "        t['S2CN1'] = t[t.index == sp[1]]['CandidatesN'].values[0]\n",
    "        t['S2CAL1'] = t[t.index == sp[1]]['CandidatesAL'].values[0]\n",
    "        t['S2CN2'] = t[t.index == sp[2]]['CandidatesN'].values[0]\n",
    "        t['S2CAL2'] = t[t.index == sp[2]]['CandidatesAL'].values[0]\n",
    "\n",
    "    ts.append(t)\n",
    "\n",
    "t = pd.concat(ts).reset_index(drop=True)\n",
    "plt.figure(figsize=(5, 4))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(data=t, x='CandidatesN', y='CandidatesAL',\n",
    "             hue='modification',\n",
    "             hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "             sort=False,\n",
    "             estimator=None,\n",
    "             alpha=0.5,\n",
    "             ax=ax)\n",
    "\n",
    "# sns.scatterplot(data=t, x='S2CN', y='S2CAL',\n",
    "#                 hue='modification',\n",
    "#                 hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "#                 alpha=0.5, size=10, legend=False,\n",
    "#                 ax=ax)\n",
    "sns.scatterplot(data=t, x='S2CN0', y='S2CAL0', color='green',\n",
    "             alpha=0.3, size=10, legend=False,\n",
    "             ax=ax)\n",
    "sns.scatterplot(data=t, x='S2CN1', y='S2CAL1', color='red',\n",
    "                alpha=0.3, size=10, legend=False,\n",
    "                ax=ax)\n",
    "sns.scatterplot(data=t, x='S2CN2', y='S2CAL2', color='blue',\n",
    "                alpha=0.3, size=10, legend=False,\n",
    "                ax=ax)\n",
    "\n",
    "ax.set(xscale='log')\n",
    "ax.set(yscale='log')\n",
    "ax.set_title('Candidates number vs length by sensitivity')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "d96f29184dfdd3c1",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sensitivities = list(sorted(full_sensitivity['Sensitivity'].unique()))\n",
    "ts = []\n",
    "for (m, c, r, ds) in EXAMPLES:\n",
    "    t = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                         (full_sensitivity['modification'] == m) &\n",
    "                         (full_sensitivity['cell'] == c) &\n",
    "                         (full_sensitivity['replicate'] == r) &\n",
    "                         (full_sensitivity['Gap'] == 0)].copy().reset_index(drop=True)\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    sp = sps[name] if name in sps else None\n",
    "    t['S2N'] = sp[1]\n",
    "    t['S2CN'] = t[t.index == sp[1]]['CandidatesN'].values[0]\n",
    "    t['S2CAL'] = t[t.index == sp[1]]['CandidatesAL'].values[0]\n",
    "    ts.append(t)\n",
    "t = pd.concat(ts).reset_index(drop=True)\n",
    "display(t.sample(10))\n"
   ],
   "id": "38d5ea6957190e43",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for y in ['CandidatesN', 'CandidatesAL', 'SignalNoiseRatio', 'SignalControlRatio']:\n",
    "    print(y)\n",
    "    plt.figure(figsize=(5, 2))\n",
    "    ax = plt.axes()\n",
    "    sns.lineplot(data=t, x='SensitivityN', y=y, hue='modification',\n",
    "                 hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "                 alpha=0.5,\n",
    "                 sort=False,\n",
    "                 estimator=None,\n",
    "                 ax=ax)\n",
    "    if y == 'CandidatesN':\n",
    "        sns.scatterplot(data=t, x='S2N', y='S2CN',\n",
    "                        hue='modification',\n",
    "                        hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "                        alpha=0.8, size=50, legend=False,\n",
    "                        ax=ax)\n",
    "    elif y == 'CandidatesAL':\n",
    "        sns.scatterplot(data=t, x='S2N', y='S2CAL',\n",
    "                        hue='modification',\n",
    "                        hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "                        alpha=0.8, size=50, legend=False,\n",
    "                        ax=ax)\n",
    "    ax.set(yscale='log')\n",
    "    ax.set_title(y)\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n"
   ],
   "id": "50bde0ac3469f6b4",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(data=t, x='SensitivityN', y='CandidatesN', hue='modification',\n",
    "             hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "             alpha=0.5,\n",
    "             sort=False,\n",
    "             estimator=None,\n",
    "             ax=ax)\n",
    "ax.set(yscale='log')\n",
    "ax.set_title('Number of candidates')\n",
    "ax.set_xlabel('PEP threshold rank')\n",
    "ax.set_ylabel('Number')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ],
   "id": "4f466dc64d0c937d",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Fragmentation",
   "id": "48c1c19c9283f4a3"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for (ds, m, c, r), t in full_sensitivity.groupby(['dataset', 'modification', 'cell', 'replicate']):\n",
    "    if len(t['file'].unique()) != 1:\n",
    "        continue\n",
    "    gaps = t['file'].values[0].replace('.sensitivity.tsv', '.gaps.tsv')\n",
    "    if not os.path.exists(gaps):\n",
    "        continue\n",
    "    print(ds, m, c, r, gaps)\n",
    "    t = pd.read_csv(gaps, sep='\\t')\n",
    "    t['dataset'] = ds\n",
    "    t['modification'] = m\n",
    "    t['cell'] = c\n",
    "    t['replicate'] = r\n",
    "    t['N'] = t['CandidatesN'] / np.max(t['CandidatesN'])\n",
    "    ts.append(t)\n",
    "gaps_df = pd.concat(ts).reset_index(drop=True)\n",
    "del ts\n",
    "gaps_df.sample(5)"
   ],
   "id": "fee0cc71d4725bd6",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## From files",
   "id": "43ecf18df198c6a7"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "FILE_EXAMPLES = [\n",
    "    ('CTCF', 'H1', ''),\n",
    "    ('ATAC', 'Blood_Monocyte', 'SRR5799491', 'ImmGen'),\n",
    "    # ('ATAC', 'Blood_Monocyte', 'SRR5799492', 'ImmGen'),\n",
    "    # ('ATAC', 'Small_Intestine_CD8_T_Cell', 'SRR5799444', 'ImmGen'),\n",
    "    # ('ATAC', 'Small_Intestine_CD8_T_Cell', 'SRR5799445', 'ImmGen'),\n",
    "    # ('H3K4me3', 'CD34', '', 'Roadmap'),\n",
    "    # ('H3K4me3', 'K562', 'rep1', 'ENCODE'),\n",
    "    # ('H3K4me3', 'CD34', 'rep2', 'Roadmap'),\n",
    "    # ('H3K4me3', 'CD4', '', 'Roadmap'),\n",
    "    # ('H3K27ac', 'K562', 'rep1', 'ENCODE'),\n",
    "    # ('H3K4me1', 'H1', 'rep1', 'ENCODE'),\n",
    "    # ('H3K4me1', 'H1', 'rep2', 'ENCODE'),\n",
    "    # ('H3K27me3', 'GM12878', 'rep2', 'ENCODE'),\n",
    "    # ('H3K27me3', 'NHLF', 'rep1', 'ENCODE'),\n",
    "    # ('H3K27me3', 'CD4', '', 'Roadmap'),\n",
    "    # ('H3K27me3', 'CD4', 'rep2', 'Roadmap'),\n",
    "    ('H3K27me3', 'Monocyte', '', 'Roadmap'),\n",
    "    # ('H3K27me3', 'BCell', '', 'Roadmap'),\n",
    "    # ('H3K27me3', 'PBMC', '', 'Roadmap'),\n",
    "    # ('H3K27me3', 'TCell', '', 'Roadmap'),\n",
    "    # ('H3K27me3', 'H1', 'rep1', 'ENCODE'),\n",
    "    # ('H3K27me3', 'HepG2', 'rep1', 'ENCODE'),\n",
    "    # ('H3K27me3', 'HepG2', 'rep2', 'ENCODE'),\n",
    "    # ('H3K36me3', 'CD4', '', 'Roadmap'),\n",
    "    # ('H3K36me3', 'CD4', 'rep1', 'Roadmap'),\n",
    "    # ('H3K36me3', 'GM12878', 'rep1', 'ENCODE'),\n",
    "    # ('H3K36me3', 'CD34', '', 'Roadmap'),\n",
    "    # ('H3K36me3', 'Monocyte', '', 'Roadmap'),\n",
    "    # ('H3K36me3', 'BCell', '', 'Roadmap'),\n",
    "    # ('H3K36me3', 'TCell', '', 'Roadmap'),\n",
    "    # ('H3K36me3', 'PBMC', '', 'Roadmap'),\n",
    "]"
   ],
   "id": "4a00a85a793180a6",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r, ds) in FILE_EXAMPLES:\n",
    "    t = gaps_df[(gaps_df['dataset'] == ds) &\n",
    "                (gaps_df['modification'] == m) &\n",
    "                (gaps_df['cell'] == c) &\n",
    "                (gaps_df['replicate'] == r)]\n",
    "    if len(t) == 0:\n",
    "        continue\n",
    "    plt.figure(figsize=(5, 2))\n",
    "    ax = plt.axes()\n",
    "    sns.lineplot(data=t.reset_index(drop=True), x='Gap', y='N',\n",
    "                 estimator=None,\n",
    "                 alpha=0.5,\n",
    "                 sort=False,\n",
    "                 ax=ax)\n",
    "    ax.set_ylabel('%')\n",
    "    ax.set_title(f'{ds} {m} {c} {r}')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "8e255b345e254c97",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Per modification",
   "id": "5eae4bcd041196ba"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(data=gaps_df, x='Gap', y='N', hue='modification',\n",
    "             alpha=0.5,\n",
    "             sort=False, errorbar=('se', .95),\n",
    "             hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "             ax=ax)\n",
    "# ax.axvline(x=20, ymin=60, ymax=100, color='black', lw=1, ls='dotted')\n",
    "# sns.scatterplot(x=[20], y=[75],\n",
    "#                 color='red', s=20,\n",
    "#                 legend=False, ax=ax)\n",
    "ax.set_ylabel('Fraction')\n",
    "# ax.set_ylim(20, 100)\n",
    "# ax.set(xscale='log')\n",
    "ax.set_title('Fragmentation')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "769ba68235053c78",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.rc('font', size=5)\n",
    "plt.figure(figsize=(3, 6))\n",
    "axs = [plt.subplot(7, 1, i + 1) for i in range(7)]\n",
    "ymin = np.min(gaps_df['N'])\n",
    "for i, m in enumerate(['CTCF', 'ATAC'] + MODIFICATIONS):\n",
    "    t = gaps_df[gaps_df['modification'] == m]\n",
    "    if len(t) == 0:\n",
    "        continue\n",
    "    ax = axs[i]\n",
    "    sns.lineplot(t,\n",
    "                 x='Gap', y='N',\n",
    "                 hue='dataset',\n",
    "                 alpha=0.5, errorbar=('se', .95),\n",
    "                 ax=ax)\n",
    "    # sns.scatterplot(x=[20], y=[70],\n",
    "    #                 color='red', s=20,\n",
    "    #                 legend=False, ax=ax)\n",
    "    ax.set_ylim(bottom=ymin)\n",
    "    ax.set_ylabel('Fraction')\n",
    "    ax.set_title(f'{m} fragmentation')\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "plt.rc('font', size=8)"
   ],
   "id": "797e5ca422a4064d",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Autocorrelations",
   "id": "2e85e862ee4771be"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Coverage",
   "id": "4621aedea8d8af90"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for m in ['CTCF', 'ATAC'] + MODIFICATIONS:\n",
    "    for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "        print(ds, m)\n",
    "        dfdsm = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                                 (full_sensitivity['modification'] == m) &\n",
    "                                 (full_sensitivity['Gap'] == 0)]\n",
    "        for (cell, rep), t in dfdsm.groupby(['cell', 'replicate']):\n",
    "            if len(t['file'].unique()) != 1:\n",
    "                continue\n",
    "            corrs = t['file'].values[0].replace('.sensitivity.tsv', '.ac.coverage.tsv')\n",
    "            if not os.path.exists(corrs):\n",
    "                continue\n",
    "            t = pd.read_csv(corrs, sep='\\t')\n",
    "            t['dataset'] = ds\n",
    "            t['modification'] = m\n",
    "            t['cell'] = cell\n",
    "            t['replicate'] = rep\n",
    "            ts.append(t)\n",
    "ac_cov_df = pd.concat(ts).reset_index(drop=True)\n",
    "del ts\n",
    "ac_cov_df.sample(5)"
   ],
   "id": "833399d4d735c379",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(ac_cov_df[ac_cov_df['D'] <= 200],\n",
    "             x='D', y='Correlation',\n",
    "             hue='modification',\n",
    "             hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "             alpha=0.5,  # errorbar=('se', .95),\n",
    "             ax=ax)\n",
    "ax.set_xlabel('Distance')\n",
    "ax.set_ylabel('Autocorrelation')\n",
    "ax.set_title('Raw signal coverage autocorrelation')\n",
    "ax.set(xscale='log')\n",
    "# ax.invert_xaxis()\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "f69fb6e48d9855ad",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.rc('font', size=5)\n",
    "plt.figure(figsize=(3, 6))\n",
    "axs = [plt.subplot(7, 1, i + 1) for i in range(7)]\n",
    "ymin = np.min(ac_cov_df['Correlation'])\n",
    "for i, m in enumerate(['CTCF', 'ATAC'] + MODIFICATIONS):\n",
    "    t = ac_cov_df[ac_cov_df['modification'] == m]\n",
    "    if len(t) == 0:\n",
    "        continue\n",
    "    ax = axs[i]\n",
    "    sns.lineplot(t,\n",
    "                 x='D', y='Correlation',\n",
    "                 hue='dataset',\n",
    "                 alpha=0.5,\n",
    "                 ax=ax)\n",
    "    ax.set_title(f'{m} Coverage autocorrelation')\n",
    "    ax.set(xscale='log')\n",
    "    ax.set_ylim(bottom=ymin)\n",
    "    # ax.invert_xaxis()\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "plt.rc('font', size=8)"
   ],
   "id": "883af62673f775aa",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## PVals",
   "id": "c4fa440bb7dbb0c3"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for m in ['CTCF', 'ATAC'] + MODIFICATIONS:\n",
    "    for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "        print(ds, m)\n",
    "        dfdsm = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                                 (full_sensitivity['modification'] == m) &\n",
    "                                 (full_sensitivity['Gap'] == 0)]\n",
    "        for (cell, rep), t in dfdsm.groupby(['cell', 'replicate']):\n",
    "            if len(t['file'].unique()) != 1:\n",
    "                continue\n",
    "            corrs = t['file'].values[0].replace('.sensitivity.tsv', '.ac.pvals.tsv')\n",
    "            if not os.path.exists(corrs):\n",
    "                continue\n",
    "            t = pd.read_csv(corrs, sep='\\t')\n",
    "            t['dataset'] = ds\n",
    "            t['modification'] = m\n",
    "            t['cell'] = cell\n",
    "            t['replicate'] = rep\n",
    "            ts.append(t)\n",
    "ac_pvals_df = pd.concat(ts).reset_index(drop=True)\n",
    "del ts\n",
    "ac_pvals_df.sample(5)"
   ],
   "id": "61a948932cdc0e9d",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(ac_pvals_df[ac_pvals_df['D'] <= 200],\n",
    "             x='D', y='Correlation',\n",
    "             hue='modification',\n",
    "             hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "             alpha=0.5,  #errorbar=('se', .95),\n",
    "             ax=ax)\n",
    "# ax.axvline(x=10, ymin=0, ymax=1.0, color='black', lw=.5, ls='dotted')\n",
    "# sns.scatterplot(x=[10], y=[0.7],\n",
    "#                 color='red', s=20,\n",
    "#                 legend=False, ax=ax)\n",
    "ax.set_xlabel('Distance')\n",
    "ax.set_ylabel('Autocorrelation')\n",
    "ax.set_title('HMM PEP autocorrelation')\n",
    "ax.set(xscale='log')\n",
    "# ax.invert_xaxis()\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "db46c16635372874",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.rc('font', size=5)\n",
    "plt.figure(figsize=(3, 6))\n",
    "axs = [plt.subplot(7, 1, i + 1) for i in range(7)]\n",
    "ymin = np.min(ac_pvals_df['Correlation'])\n",
    "for i, m in enumerate(['CTCF', 'ATAC'] + MODIFICATIONS):\n",
    "    t = ac_pvals_df[ac_pvals_df['modification'] == m]\n",
    "    if len(t) == 0:\n",
    "        continue\n",
    "    ax = axs[i]\n",
    "    sns.lineplot(t,\n",
    "                 x='D', y='Correlation',\n",
    "                 hue='dataset',\n",
    "                 alpha=0.5,  #errorbar=('se', .95),\n",
    "                 ax=ax)\n",
    "    # sns.scatterplot(x=[10], y=[0.7],\n",
    "    #                 color='red', s=20,\n",
    "    #                 legend=False, ax=ax)\n",
    "    ax.set_ylim(bottom=ymin)\n",
    "    ax.set_title(f'{m} Pvals autocorrelation')\n",
    "    ax.set(xscale='log')\n",
    "    # ax.invert_xaxis()\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "plt.rc('font', size=8)"
   ],
   "id": "abbf209487e819ee",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Per modification",
   "id": "f9a5ea6eeba6d224"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# import plotly.express as px\n",
    "#\n",
    "# ac_pvals_df['experiment'] = ac_pvals_df['dataset'] + ' ' + ac_pvals_df['cell'] + ' ' + ac_pvals_df['replicate']\n",
    "#\n",
    "# plt.rc('font', size=4)\n",
    "# # for m in ['ATAC'] + MODIFICATIONS:\n",
    "# for m in ['H3K27me3']:\n",
    "#     dfm = ac_pvals_df[(ac_pvals_df['modification'] == m) & (ac_pvals_df['D'] <= 50)]\n",
    "#     if len(dfm) == 0:\n",
    "#         continue\n",
    "#     fig = px.line(dfm, x='D', y='Correlation', color='experiment',\n",
    "#                   hover_data=['experiment'],\n",
    "#                   title=f'{m} Pvals autocorrelation',\n",
    "#                   width=1000, height=400)\n",
    "#     fig.show()\n",
    "# plt.rc('font', size=8)"
   ],
   "id": "9d8685e5f8e24edd",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Candidates",
   "id": "927535d8f0a82efc"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "cns = {}\n",
    "cals = {}\n",
    "for (ds, m, c, r), t in full_sensitivity.groupby(['dataset', 'modification', 'cell', 'replicate']):\n",
    "    tt = t[t['Gap'] == 0].sort_values('Sensitivity').reset_index(drop=True)\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    sp = sps[name]\n",
    "    cn = tt['CandidatesN'].values[sp[1]]\n",
    "    cal = tt['CandidatesAL'].values[sp[1]]\n",
    "    cns[(ds, m, c, r)] = cn\n",
    "    cals[(ds, m, c, r)] = cal\n",
    "\n",
    "full_infos['candidatesN'] = [cns[(ds, m, c, r)] if (ds, m, c, r) in cns else 0 for _, (ds, m, c, r) in\n",
    "                             full_infos[['dataset', 'modification', 'cell', 'replicate']].iterrows()]\n",
    "full_infos['candidatesAL'] = [cals[(ds, m, c, r)] if (ds, m, c, r) in cals else 0 for _, (ds, m, c, r) in\n",
    "                              full_infos[['dataset', 'modification', 'cell', 'replicate']].iterrows()]\n",
    "full_infos.sample(3)"
   ],
   "id": "ba9ab4a704259a17",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Remove outliers\n",
    "OUTLIERS = [\n",
    "    ('H3K4me3', 'NK', ''),\n",
    "    ('H3K27me3', 'TCell', ''),\n",
    "    ('H3K27me3', 'BCell', ''),\n",
    "]\n",
    "df_peaks = full_infos.loc[[(m, c, r) not in OUTLIERS\n",
    "                           for _, (m, c, r) in full_infos[['modification', 'cell', 'replicate']].iterrows()]].copy()"
   ],
   "id": "9adef09000c4c46b",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sns.set_style('whitegrid')\n",
    "plt.figure(figsize=(5, 3))\n",
    "ax = plt.axes()\n",
    "t = df_peaks\n",
    "sns.barplot(data=t, x='modification', y='candidatesN', hue='dataset',\n",
    "            capsize=.1, edgecolor=\"black\",\n",
    "            order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "            hue_order=['ENCODE', 'Roadmap', 'ABF', 'ImmGen'],\n",
    "            err_kws={'linewidth': 2},\n",
    "            ax=ax)\n",
    "sns.stripplot(data=t, x='modification', y='candidatesN', hue='dataset',\n",
    "              dodge=True, size=2, palette='dark:black', alpha=0.5,\n",
    "              order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "              hue_order=['ENCODE', 'Roadmap', 'ABF', 'ImmGen'], legend=False,\n",
    "              ax=ax)\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "# ax.set_ylim(0, 2e5)\n",
    "# ax.set(yscale='log')\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()"
   ],
   "id": "3bb5baa8b8cc94b2",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "full_infos[(full_infos['modification'] == 'H3K27me3')].sort_values(by=['candidatesN'], ascending=False)[\n",
    "    ['cell', 'replicate', 'candidatesN']].head(10)"
   ],
   "id": "bba47f0f05495e53",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Min signal-to-noise",
   "id": "4a0a10d347806057"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "min_sns = {}\n",
    "for m in ['CTCF', 'ATAC'] + MODIFICATIONS:\n",
    "    for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "        print(ds, m)\n",
    "        dfdsm = full_sensitivity[(full_sensitivity['dataset'] == ds) & (full_sensitivity['modification'] == m)]\n",
    "        for (cell, rep), t in dfdsm.groupby(['cell', 'replicate']):\n",
    "            tt = t[t['Gap'] == 0].sort_values(by=['Sensitivity']).copy().reset_index(drop=True)\n",
    "            sp = detect_sensitivity_triangle(tt)\n",
    "            snr = tt['SignalNoiseRatio'].copy().values\n",
    "            snr[:sp[0]] = 1e6\n",
    "            snr[sp[2]:] = 1e6\n",
    "            mi = np.argmin(snr)\n",
    "            min_snr = snr[mi]\n",
    "            print(f'{ds} {m} {cell} {rep} {mi}: {min_snr}')\n",
    "            min_sns[(ds, m, cell, rep)] = min_snr\n",
    "\n",
    "full_infos['minimal_signal_noise_ratio'] = [min_sns[(ds, m, c, r)] for _, (ds, m, c, r) in\n",
    "                                            full_infos[['dataset', 'modification', 'cell', 'replicate']].iterrows()]\n",
    "full_infos.sample(3)"
   ],
   "id": "5810bd0c815a4042",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Pvalues",
   "id": "cd31dab11b79b5a2"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Examples",
   "id": "6ac6727a728149ef"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "t = StringIO('''Q\tLogNullP\n",
    "0.0\t-216.5814666748047\n",
    "1.0E-5\t-119.61390686035156\n",
    "2.0E-5\t-112.77011108398438\n",
    "3.0000000000000004E-5\t-105.98941040039062\n",
    "4.0E-5\t-102.6241455078125\n",
    "5.0E-5\t-102.6241455078125\n",
    "6.0E-5\t-99.27644348144531\n",
    "7.000000000000001E-5\t-99.27644348144531\n",
    "8.0E-5\t-95.94697570800781\n",
    "9.0E-5\t-92.6364517211914\n",
    "1.0E-4\t-92.6364517211914\n",
    "1.1E-4\t-92.6364517211914\n",
    "1.2E-4\t-89.34564208984375\n",
    "1.3000000000000002E-4\t-89.34564208984375\n",
    "1.4000000000000001E-4\t-89.34564208984375\n",
    "1.5000000000000001E-4\t-86.07534790039062\n",
    "1.6E-4\t-86.07534790039062\n",
    "1.7E-4\t-86.07534790039062\n",
    "1.8E-4\t-86.07534790039062\n",
    "1.9E-4\t-86.07534790039062\n",
    "2.0E-4\t-82.82644653320312\n",
    "2.1E-4\t-82.82644653320312\n",
    "2.2E-4\t-82.82644653320312\n",
    "2.3E-4\t-82.82644653320312\n",
    "2.4E-4\t-82.82644653320312\n",
    "2.5E-4\t-79.59986877441406\n",
    "2.6000000000000003E-4\t-79.59986877441406\n",
    "2.7000000000000006E-4\t-79.59986877441406\n",
    "2.800000000000001E-4\t-79.59986877441406\n",
    "2.900000000000001E-4\t-79.59986877441406\n",
    "3.0000000000000014E-4\t-79.59986877441406\n",
    "3.1000000000000016E-4\t-76.39661407470703\n",
    "3.200000000000002E-4\t-76.39661407470703\n",
    "3.300000000000002E-4\t-76.39661407470703\n",
    "3.4000000000000024E-4\t-76.39661407470703\n",
    "3.5000000000000027E-4\t-76.39661407470703\n",
    "3.600000000000003E-4\t-76.39661407470703\n",
    "3.700000000000003E-4\t-76.39661407470703\n",
    "3.8000000000000035E-4\t-76.3956298828125\n",
    "3.9000000000000037E-4\t-73.21774291992188\n",
    "4.000000000000004E-4\t-73.21774291992188\n",
    "4.100000000000004E-4\t-73.21774291992188\n",
    "4.2000000000000045E-4\t-73.21774291992188\n",
    "4.300000000000005E-4\t-73.21774291992188\n",
    "4.400000000000005E-4\t-73.21774291992188\n",
    "4.5000000000000053E-4\t-73.21774291992188\n",
    "4.6000000000000056E-4\t-73.21774291992188\n",
    "4.700000000000006E-4\t-70.0644302368164\n",
    "4.800000000000006E-4\t-70.0644302368164\n",
    "4.900000000000006E-4\t-70.0644302368164\n",
    "5.000000000000007E-4\t-70.0644302368164\n",
    "5.100000000000007E-4\t-70.0644302368164\n",
    "5.200000000000007E-4\t-70.0644302368164\n",
    "5.300000000000007E-4\t-70.0644302368164\n",
    "5.400000000000008E-4\t-70.0644302368164\n",
    "5.500000000000008E-4\t-70.0644302368164\n",
    "5.600000000000008E-4\t-70.0644302368164\n",
    "5.700000000000008E-4\t-70.0644302368164\n",
    "5.800000000000009E-4\t-68.97872924804688\n",
    "5.900000000000009E-4\t-66.93791198730469\n",
    "6.000000000000009E-4\t-66.93791198730469\n",
    "6.10000000000001E-4\t-66.93791198730469\n",
    "6.20000000000001E-4\t-66.93791198730469\n",
    "6.30000000000001E-4\t-66.93791198730469\n",
    "6.40000000000001E-4\t-66.93791198730469\n",
    "6.500000000000011E-4\t-66.93791198730469\n",
    "6.600000000000011E-4\t-66.93791198730469\n",
    "6.700000000000011E-4\t-66.93791198730469\n",
    "6.800000000000011E-4\t-66.93791198730469\n",
    "6.900000000000012E-4\t-66.93791198730469\n",
    "7.000000000000012E-4\t-66.93778228759766\n",
    "7.100000000000012E-4\t-63.83954620361328\n",
    "7.200000000000012E-4\t-63.83954620361328\n",
    "7.300000000000013E-4\t-63.83954620361328\n",
    "7.400000000000013E-4\t-63.83954620361328\n",
    "7.500000000000013E-4\t-63.83954620361328\n",
    "7.600000000000013E-4\t-63.83954620361328\n",
    "7.700000000000014E-4\t-63.83954620361328\n",
    "7.800000000000014E-4\t-63.83954620361328\n",
    "7.900000000000014E-4\t-63.83954620361328\n",
    "8.000000000000014E-4\t-63.83954620361328\n",
    "8.100000000000015E-4\t-63.83954620361328\n",
    "8.200000000000015E-4\t-63.83954620361328\n",
    "8.300000000000015E-4\t-63.83954620361328\n",
    "8.400000000000016E-4\t-63.83954620361328\n",
    "8.500000000000016E-4\t-63.8370475769043\n",
    "8.600000000000016E-4\t-60.77081298828125\n",
    "8.700000000000016E-4\t-60.77081298828125\n",
    "8.800000000000017E-4\t-60.77081298828125\n",
    "8.900000000000017E-4\t-60.77081298828125\n",
    "9.000000000000017E-4\t-60.77081298828125\n",
    "9.100000000000017E-4\t-60.77081298828125\n",
    "9.200000000000018E-4\t-60.77081298828125\n",
    "9.300000000000018E-4\t-60.77081298828125\n",
    "9.400000000000018E-4\t-60.77081298828125\n",
    "9.500000000000018E-4\t-60.77081298828125\n",
    "9.600000000000019E-4\t-60.77081298828125\n",
    "9.700000000000019E-4\t-60.77081298828125\n",
    "9.80000000000002E-4\t-60.77081298828125\n",
    "9.90000000000002E-4\t-60.77081298828125\n",
    "''')\n",
    "t = pd.read_csv(t, sep='\\\\t')\n",
    "plt.figure(figsize=(5, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(data=t, x='Q', y='LogNullP', estimator=None,\n",
    "             alpha=0.5,\n",
    "             sort=False,\n",
    "             ax=ax)\n",
    "ax.set_title('K27me3 Monocyte')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "22426ecde2e4bc7d",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## From files",
   "id": "a4c80049b250fd0d"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for m in ['CTCF', 'ATAC'] + MODIFICATIONS:\n",
    "    for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "        print(ds, m)\n",
    "        dfdsm = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                                 (full_sensitivity['modification'] == m) &\n",
    "                                 (full_sensitivity['Gap'] == 0)]\n",
    "        for (cell, rep), t in dfdsm.groupby(['cell', 'replicate']):\n",
    "            if len(t['file'].unique()) != 1:\n",
    "                continue\n",
    "            pvals = t['file'].values[0].replace('.sensitivity.tsv', '.logps.tsv')\n",
    "            if not os.path.exists(pvals):\n",
    "                continue\n",
    "            t = pd.read_csv(pvals, sep='\\t')\n",
    "            t['dataset'] = ds\n",
    "            t['modification'] = m\n",
    "            t['cell'] = cell\n",
    "            t['replicate'] = rep\n",
    "            ts.append(t)\n",
    "pvals_df = pd.concat(ts).reset_index(drop=True)\n",
    "del ts\n",
    "pvals_df.sample(5)"
   ],
   "id": "34bba4673343c4e5",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r, ds) in FILE_EXAMPLES:\n",
    "    t = pvals_df[(pvals_df['dataset'] == ds) &\n",
    "                 (pvals_df['modification'] == m) &\n",
    "                 (pvals_df['cell'] == c) &\n",
    "                 (pvals_df['replicate'] == r)].copy().reset_index(drop=True)\n",
    "    if len(t) == 0:\n",
    "        continue\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    plt.figure(figsize=(4, 2))\n",
    "    ax = plt.axes()\n",
    "    sns.lineplot(t,\n",
    "                 x='Q', y='LogNullP', estimator=None,\n",
    "                 ax=ax)\n",
    "    ax.set_title(f'Pvalues {name}')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "c7df342e6742d71e",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## All",
   "id": "4182e9bc040be1bc"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "pvals_df['MinusLogNullP'] = -pvals_df['LogNullP']\n",
    "plt.figure(figsize=(5, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(pvals_df[pvals_df['Q']<=0.005],\n",
    "             x='Q', y='LogNullP',\n",
    "             hue='modification',\n",
    "             hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "             alpha=0.5, errorbar=('se', .95),\n",
    "             ax=ax)\n",
    "ax.set_title('Pvalues')\n",
    "# ax.set(yscale='log')\n",
    "# ax.invert_xaxis()\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "9c36a2f611cc35b2",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.rc('font', size=5)\n",
    "plt.figure(figsize=(3, 6))\n",
    "axs = [plt.subplot(7, 1, i + 1) for i in range(7)]\n",
    "for i, m in enumerate(['CTCF', 'ATAC'] + MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    sns.lineplot(pvals_df[(pvals_df['modification'] == m) & (pvals_df['dataset'] != '') & (pvals_df['Q'] <= 0.005)],\n",
    "                 x='Q', y='LogNullP',\n",
    "                 hue='dataset',\n",
    "                 alpha=0.5, errorbar=('se', .95),\n",
    "                 ax=ax)\n",
    "    ax.set_title(f'{m} Pvalues')\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "plt.rc('font', size=8)"
   ],
   "id": "b36966de8e396dc4",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "5d849a66e2dd6229",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Summary",
   "id": "47421e63e17b67b"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "full_infos['effective_coverage'] = full_infos['coverage'] - full_infos['beta'] * full_infos['control_coverage'] * \\\n",
    "                                   full_infos['control_scale']"
   ],
   "id": "fd7bd1fd3fff66e6",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "features = [\n",
    "    'model_signal',\n",
    "    'model_noise',\n",
    "    'model_signal_to_noise',\n",
    "    'log_null_pvals_mean',\n",
    "    'log_null_pvals_std',\n",
    "    'coverage',\n",
    "    'control_coverage',\n",
    "    'control_scale',\n",
    "    'beta',\n",
    "    'min_correlation',\n",
    "    'coverage_non_zero',\n",
    "    'coverage_max',\n",
    "    'coverage_mean',\n",
    "    'coverage_median',\n",
    "    'coverage_std',\n",
    "    # 'roughness',\n",
    "    'signal_density',\n",
    "    'noise_density',\n",
    "    'signal_to_noise',\n",
    "    'signal_to_control',\n",
    "    'autocorrelation_average_score',\n",
    "    'fragmentation_average_score'\n",
    "]\n",
    "\n",
    "full_features = features # + ['candidatesN', 'candidatesAL', 'minimal_signal_noise_ratio', ]"
   ],
   "id": "2d90931a0b74fa89",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "full_infos[(full_infos['modification'] == 'H3K27ac') &\n",
    "           (full_infos['cell'].isin(['CD4']))][['modification', 'cell'] + full_features].T"
   ],
   "id": "75b9688bb3c67953",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ncols = int(math.ceil(math.sqrt(len(full_features) + 3)))\n",
    "nrows = int(math.ceil(len(full_features) / ncols))\n",
    "plt.figure(figsize=(3.2 * ncols, 2.5 * nrows))\n",
    "axs = [plt.subplot(nrows, ncols, i + 1) for i in range(len(full_features))]\n",
    "for i, y in enumerate(full_features):\n",
    "    ax = axs[i]\n",
    "    sns.boxplot(data=full_infos, x='modification', y=y, \n",
    "                hue='dataset',\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF', 'ImmGen'],\n",
    "                order=['ATAC'] + MODIFICATIONS,\n",
    "                showfliers=False,\n",
    "                ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\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()"
   ],
   "id": "53299009085c235a",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "full_infos.fillna(0, inplace=True)\n",
    "corrs = full_infos[full_features].corr(method='pearson').fillna(0)\n",
    "plt.rc('font', size=5)\n",
    "sns.clustermap(corrs,\n",
    "               col_cluster=True, row_cluster=True,\n",
    "               figsize=(4.5, 4.5),\n",
    "               cmap=plt.cm.seismic,\n",
    "               vmin=-1, vmax=1)\n",
    "plt.show()\n",
    "plt.rc('font', size=8)"
   ],
   "id": "24c2ed734ad1a421",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# for m in MODIFICATIONS:\n",
    "#     print(m)\n",
    "#     corrs = full_infos[full_infos['modification'] == m][features].corr(method='spearman')\n",
    "#     plt.rc('font', size=5)\n",
    "#     sns.clustermap(corrs,\n",
    "#                    col_cluster=True, row_cluster=True,\n",
    "#                    figsize=(4, 4),\n",
    "#                    cmap=plt.cm.seismic,\n",
    "#                 vmin=-1, vmax=1)\n",
    "#     plt.show()\n",
    "#     plt.rc('font', size=8)"
   ],
   "id": "216a45728b85adf3",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Segments",
   "id": "2b98ee6a52fadda1"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import matplotlib.patches as mpatches\n",
    "\n",
    "\n",
    "def d(a, b):\n",
    "    return a / b if b != 0 else 0\n",
    "\n",
    "\n",
    "def plot_segments(sdf, df, name, sp=None, ax=None, legend=True):\n",
    "    sdf['New %'] = [d(n, t) * 100 for n, t in zip(sdf['New'], sdf['CandidatesN'])]\n",
    "    sdf['Total %'] = 100\n",
    "    axisnone = ax is None\n",
    "    if ax is None:\n",
    "        plt.figure(figsize=(8, 2))\n",
    "        ax = plt.axes()\n",
    "    # Background\n",
    "    sns.barplot(x=range(len(sdf)), y=sdf['Total %'],\n",
    "                color='orange',\n",
    "                alpha=0.5,\n",
    "                ax=ax)\n",
    "    # Foreground\n",
    "    sns.barplot(x=range(len(sdf)), y=sdf['New %'],\n",
    "                color='blue',\n",
    "                alpha=0.5,\n",
    "                ax=ax)\n",
    "    if df is not None:\n",
    "        assert len(sdf) == len(df)\n",
    "        df = df.sort_values(by=['Sensitivity'])\n",
    "        ns = np.log1p(df['CandidatesN'])\n",
    "        ns = ns * 80 / ns.max() + 10\n",
    "        sns.lineplot(x=range(len(sdf)), y=ns,\n",
    "                     color='black',\n",
    "                     sort=False,\n",
    "                     estimator=None,\n",
    "                     ax=ax)\n",
    "        als = np.log1p(df['CandidatesAL'])\n",
    "        als = als * 80 / als.max() + 10\n",
    "        sns.lineplot(x=range(len(sdf)), y=als,\n",
    "                     color='brown',\n",
    "                     sort=False,\n",
    "                     estimator=None,\n",
    "                     ax=ax)\n",
    "        # snrs = np.log1p(df['SignalNoiseRatio'])\n",
    "        # snrs = snrs * 80 / snrs.max() + 10\n",
    "        # sns.lineplot(x=range(len(sdf)), y=snrs,\n",
    "        #              color='green',\n",
    "        #              sort=False,\n",
    "        #              estimator=None,\n",
    "        #              ax=ax)\n",
    "\n",
    "    ax.set_title(name)\n",
    "    ax.set_ylabel('peaks %')\n",
    "\n",
    "    # Limit ticks number\n",
    "    # for i, label in enumerate(ax.get_xticklabels()):\n",
    "    #     if axisnone and i not in tens_closest:\n",
    "    #         label.set_visible(False)\n",
    "    sens = np.array(sorted(sdf['Sensitivity']))\n",
    "    tens_closest = {np.argmin(np.fabs(sens - x)): x for x in [-100, -10, -1, -0.1, -1e-2, -1e-3, -1e-4, -1e-6]}\n",
    "    labels = [str(tens_closest[i]) if i in tens_closest else '' for i in range(len(ax.get_xticklabels()))]\n",
    "    ax.set_xticklabels(labels)\n",
    "    ax.set_xlabel('Sensitivity')\n",
    "    for i in tens_closest.keys():\n",
    "        ax.axvline(x=i, ymin=0, ymax=90, color='grey', lw=0.5)\n",
    "\n",
    "    # Points\n",
    "    if sp is not None:\n",
    "        ax.axvline(x=sp[2], ymin=0, ymax=90, color='blue', lw=2)\n",
    "        ax.axvline(x=sp[1], ymin=0, ymax=90, color='red', lw=2)\n",
    "        ax.axvline(x=sp[0], ymin=0, ymax=90, color='green', lw=2)\n",
    "        # ax.axvline(x=int((sp[1] + sp[2]) / 2), ymin=0, ymax=90, color='grey', lw=2)\n",
    "\n",
    "    # add legend\n",
    "    if legend:\n",
    "        top_bar = mpatches.Patch(color='orange', label='old %')\n",
    "        bottom_bar = mpatches.Patch(color='blue', label='new %')\n",
    "        ns = mpatches.Patch(color='black', label='number')\n",
    "        if df is None:\n",
    "            plt.legend(handles=[bottom_bar, top_bar, ns])\n",
    "        else:\n",
    "            als = mpatches.Patch(color='brown', label='avglen')\n",
    "            # snrs = mpatches.Patch(color='green', label='snr')\n",
    "            # plt.legend(handles=[bottom_bar, top_bar, ns, als, snrs])\n",
    "            plt.legend(handles=[bottom_bar, top_bar, ns, als])\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    new_value = 1\n",
    "    for patch in ax.patches:\n",
    "        current_width = patch.get_width()\n",
    "        diff = current_width - new_value\n",
    "        # we change the bar width\n",
    "        patch.set_width(new_value)\n",
    "        # we recenter the bar\n",
    "        patch.set_x(patch.get_x() + diff * .5)\n",
    "    if axisnone:\n",
    "        plt.tight_layout()\n",
    "        plt.show()"
   ],
   "id": "244ce798ed93441c",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## From files",
   "id": "bb3845fd81690272"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r, ds) in FILE_EXAMPLES:\n",
    "    t = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                         (full_sensitivity['modification'] == m) &\n",
    "                         (full_sensitivity['cell'] == c) &\n",
    "                         (full_sensitivity['replicate'] == r) &\n",
    "                         (full_sensitivity['Gap'] == 0)]\n",
    "    if len(t['file'].unique()) != 1:\n",
    "        continue\n",
    "    segments = t['file'].values[0].replace('.sensitivity.tsv', '.segments.tsv')\n",
    "    if not os.path.exists(segments):\n",
    "        continue\n",
    "    sdf = pd.read_csv(segments, sep='\\t')\n",
    "    # Hack with indexing\n",
    "    # sdf['Sensitivity'] = sorted(sdf['Sensitivity'])\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    sp = sps[name] if name in sps else None\n",
    "    plot_segments(sdf, t, name, sp)"
   ],
   "id": "8be856cca22646ec",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r, ds) in [\n",
    "    ('H3K27me3', 'Monocyte', '', 'Roadmap'),\n",
    "]:\n",
    "    t = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                         (full_sensitivity['modification'] == m) &\n",
    "                         (full_sensitivity['cell'] == c) &\n",
    "                         (full_sensitivity['replicate'] == r) &\n",
    "                         (full_sensitivity['Gap'] == 0)].copy().reset_index(drop=True)\n",
    "    segments = t['file'].values[0].replace('.sensitivity.tsv', '.segments.tsv')\n",
    "    if not os.path.exists(segments):\n",
    "        continue\n",
    "    sdf = pd.read_csv(segments, sep='\\t')\n",
    "    sdf['New %'] = [d(n, t) * 100 for n, t in zip(sdf['New'], sdf['CandidatesN'])]\n",
    "    sdf['SensitivityN'] = range(len(sdf))\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    plt.figure(figsize=(5, 2))\n",
    "    ax = plt.axes()\n",
    "    sns.lineplot(data=sdf, x='SensitivityN', y='New %',\n",
    "             sort=False,\n",
    "             estimator=None,\n",
    "             ax=ax)\n",
    "    sens = np.array(t['Sensitivity'])\n",
    "    # First and last labels are out of the plot\n",
    "    sens_ticks = np.array([sens[int(t)] if 0 <= t < len(sens) else np.nan for t in ax.get_xticks()[1:-1]])\n",
    "    sens_ticks[-1] = sens.max()\n",
    "    labels = [''] + [f'{t:.1e}' for t in sens_ticks] + ['']\n",
    "    ax.set_xticklabels(labels)\n",
    "    ax.set_title('Percentage of new candidates')\n",
    "    ax.set_xlabel('PEP threshold')\n",
    "    ax.set_ylabel('New %    ')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n"
   ],
   "id": "cb405f23b1fe2116",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "\n",
    "for (m, c, r, ds) in EXAMPLES:\n",
    "    t = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                         (full_sensitivity['modification'] == m) &\n",
    "                         (full_sensitivity['cell'] == c) &\n",
    "                         (full_sensitivity['replicate'] == r) &\n",
    "                         (full_sensitivity['Gap'] == 0)]\n",
    "    if len(t['file'].unique()) != 1:\n",
    "        continue\n",
    "    segments = t['file'].values[0].replace('.sensitivity.tsv', '.segments.tsv')\n",
    "    if not os.path.exists(segments):\n",
    "        continue\n",
    "    sdf = pd.read_csv(segments, sep='\\t')\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    sp = sps[name] if name in sps else None\n",
    "    plot_segments(sdf, t, name, sp)"
   ],
   "id": "a9440f46f5469f5e",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Per modification",
   "id": "edee090677b35701"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for m in ['CTCF', 'ATAC'] + MODIFICATIONS:\n",
    "    for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "        print(ds, m)\n",
    "        dfdsm = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                                 (full_sensitivity['modification'] == m) &\n",
    "                                 (full_sensitivity['Gap'] == 0)]\n",
    "        for (cell, rep), t in dfdsm.groupby(['cell', 'replicate']):\n",
    "            if len(t['file'].unique()) != 1:\n",
    "                continue\n",
    "            segments = t['file'].values[0].replace('.sensitivity.tsv', '.segments.tsv')\n",
    "            if not os.path.exists(segments):\n",
    "                continue\n",
    "            sdf = pd.read_csv(segments, sep='\\t')\n",
    "            sdf['dataset'] = ds\n",
    "            sdf['modification'] = m\n",
    "            sdf['cell'] = cell\n",
    "            sdf['replicate'] = rep\n",
    "            sdf['SensitivityN'] = range(len(sdf))\n",
    "            ts.append(sdf)\n",
    "segments_df = pd.concat(ts).reset_index(drop=True)\n",
    "segments_df['New %'] = [d(n, t) * 100 for n, t in zip(segments_df['New'], segments_df['CandidatesN'])]\n",
    "del ts\n",
    "segments_df.fillna(0, inplace=True)\n",
    "segments_df.sample(5)"
   ],
   "id": "dd0d9676612a52f9",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(segments_df,\n",
    "             x='SensitivityN', y='New %',\n",
    "             hue='modification',\n",
    "             hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "             alpha=0.5, errorbar=('se', .95),\n",
    "             ax=ax)\n",
    "ax.set_xlabel('PEP threshold')\n",
    "ax.set_ylabel('New %')\n",
    "ax.set_title('Percentage of new candidates')\n",
    "# ax.set(xscale='log')\n",
    "# ax.invert_xaxis()\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "ca428de8cc452917",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.rc('font', size=5)\n",
    "plt.figure(figsize=(3, 6))\n",
    "axs = [plt.subplot(7, 1, i + 1) for i in range(7)]\n",
    "for i, m in enumerate(['CTCF', 'ATAC'] + MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    sns.lineplot(segments_df[segments_df['modification'] == m],\n",
    "                 x='SensitivityN', y='New %',\n",
    "                 hue='dataset',\n",
    "                 alpha=0.5, errorbar=('se', .95),\n",
    "                 ax=ax)\n",
    "    ax.set_ylabel('peaks %')\n",
    "    ax.set_title(f'{m} New %')\n",
    "    # ax.set(xscale='log')\n",
    "    # ax.invert_xaxis()\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "plt.rc('font', size=8)"
   ],
   "id": "1e77a5cbd0c69d79",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Individual",
   "id": "4803c50f488d9c7f"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# import plotly.express as px\n",
    "#\n",
    "# segments_df['experiment'] = segments_df['dataset'] + ' ' + segments_df['cell'] + ' ' + segments_df['replicate']\n",
    "#\n",
    "# plt.rc('font', size=4)\n",
    "# # for m in ['ATAC'] + MODIFICATIONS:\n",
    "# for m in ['ATAC']:\n",
    "#     dfm = segments_df[segments_df['modification'] == m]\n",
    "#     if len(dfm) == 0:\n",
    "#         continue\n",
    "#     fig = px.line(dfm, x='SensitivityN', y='New %', color='experiment',\n",
    "#                   hover_data=['experiment', 'Sensitivity', 'New %'],\n",
    "#                   log_x=True,\n",
    "#                   title=f'{m} New %',\n",
    "#                   width=1000, height=400)\n",
    "#     fig.show()\n",
    "# plt.rc('font', size=8)"
   ],
   "id": "e48952ae3bc4efc",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## All",
   "id": "39254ffe77c954f2"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for m in ['CTCF', 'ATAC'] + MODIFICATIONS:\n",
    "    for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "        print(ds, m)\n",
    "        dfdsm = segments_df[(segments_df['dataset'] == ds) &\n",
    "                            (segments_df['modification'] == m)]\n",
    "        n = min(len(dfdsm[['cell', 'replicate']].copy().drop_duplicates()), 18)\n",
    "        print(n)\n",
    "        if n <= 6:\n",
    "            ncols, nrows = n, 1\n",
    "        else:\n",
    "            ncols = int(math.floor(math.sqrt(n) + 2))\n",
    "            nrows = int(math.ceil(n / ncols))\n",
    "        plt.figure(figsize=(5 * ncols + 1, 2 * nrows))\n",
    "        axs = [plt.subplot(nrows, ncols, i + 1) for i in range(n)]\n",
    "        for i, ((cell, rep), sdf) in enumerate(dfdsm.groupby(['cell', 'replicate'])):\n",
    "            if i == 18:\n",
    "                break\n",
    "            t = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                                 (full_sensitivity['modification'] == m) &\n",
    "                                 (full_sensitivity['cell'] == cell) &\n",
    "                                 (full_sensitivity['replicate'] == rep) &\n",
    "                                 (full_sensitivity['Gap'] == 0)]\n",
    "            ax = axs[i]\n",
    "            name = f'{ds} {m} {cell} {rep}'\n",
    "            sp = sps[name] if name in sps else None\n",
    "            if i == n - 1:\n",
    "                plot_segments(sdf, t, name, sp, ax, True)\n",
    "                sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "            else:\n",
    "                plot_segments(sdf, t, name, sp, ax, False)\n",
    "        plt.tight_layout()\n",
    "        plt.show()"
   ],
   "id": "403b81d1b1183e46",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "638a8ae655c577e",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Coverage",
   "id": "b48a36da1525eec6"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## TCell H3K27me3 scores",
   "id": "37f65c26281e3f22"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "out = StringIO('''Q\tCoverage\n",
    "1\t1.101189395318887\n",
    "2\t1.101189395318887\n",
    "3\t1.171189395318887\n",
    "4\t1.171189395318887\n",
    "5\t1.171189395318887\n",
    "6\t1.171189395318887\n",
    "7\t1.171189395318887\n",
    "8\t1.2411893953188868\n",
    "9\t1.2411893953188868\n",
    "10\t1.2411893953188868\n",
    "11\t1.2411893953188868\n",
    "12\t1.2411893953188868\n",
    "13\t1.2411893953188868\n",
    "14\t1.2411893953188868\n",
    "15\t1.2411893953188868\n",
    "16\t1.2411893953188868\n",
    "17\t1.2411893953188868\n",
    "18\t1.2411893953188868\n",
    "19\t1.311189395318887\n",
    "20\t1.311189395318887\n",
    "21\t1.311189395318887\n",
    "22\t1.311189395318887\n",
    "23\t1.311189395318887\n",
    "24\t1.311189395318887\n",
    "25\t1.311189395318887\n",
    "26\t1.311189395318887\n",
    "27\t1.311189395318887\n",
    "28\t1.311189395318887\n",
    "29\t1.311189395318887\n",
    "30\t1.311189395318887\n",
    "31\t1.311189395318887\n",
    "32\t1.311189395318887\n",
    "33\t1.311189395318887\n",
    "34\t1.311189395318887\n",
    "35\t1.381189395318887\n",
    "36\t1.381189395318887\n",
    "37\t1.381189395318887\n",
    "38\t1.381189395318887\n",
    "39\t1.381189395318887\n",
    "40\t1.381189395318887\n",
    "41\t1.381189395318887\n",
    "42\t1.381189395318887\n",
    "43\t1.381189395318887\n",
    "44\t1.381189395318887\n",
    "45\t1.381189395318887\n",
    "46\t1.381189395318887\n",
    "47\t1.381189395318887\n",
    "48\t1.381189395318887\n",
    "49\t1.381189395318887\n",
    "50\t2.4823787906377737\n",
    "51\t2.552378790637774\n",
    "52\t2.552378790637774\n",
    "53\t2.622378790637774\n",
    "54\t2.622378790637774\n",
    "55\t2.622378790637774\n",
    "56\t2.622378790637774\n",
    "57\t2.692378790637774\n",
    "58\t2.692378790637774\n",
    "59\t2.692378790637774\n",
    "60\t2.692378790637774\n",
    "61\t2.692378790637774\n",
    "62\t2.692378790637774\n",
    "63\t2.692378790637774\n",
    "64\t2.762378790637774\n",
    "65\t2.762378790637774\n",
    "66\t2.762378790637774\n",
    "67\t2.762378790637774\n",
    "68\t2.762378790637774\n",
    "69\t2.762378790637774\n",
    "70\t2.762378790637774\n",
    "71\t3.933568185956661\n",
    "72\t4.003568185956661\n",
    "73\t4.003568185956661\n",
    "74\t4.003568185956661\n",
    "75\t4.073568185956661\n",
    "76\t4.073568185956661\n",
    "77\t4.073568185956661\n",
    "78\t4.073568185956661\n",
    "79\t4.143568185956661\n",
    "80\t4.143568185956661\n",
    "81\t4.143568185956661\n",
    "82\t5.244757581275548\n",
    "83\t5.384757581275548\n",
    "84\t5.384757581275548\n",
    "85\t5.4547575812755476\n",
    "86\t5.4547575812755476\n",
    "87\t5.524757581275548\n",
    "88\t5.524757581275548\n",
    "89\t6.695946976594435\n",
    "90\t6.8359469765944345\n",
    "91\t6.8359469765944345\n",
    "92\t6.905946976594435\n",
    "93\t8.147136371913321\n",
    "94\t8.217136371913321\n",
    "95\t8.287136371913322\n",
    "96\t9.598325767232208\n",
    "97\t10.909515162551095\n",
    "98\t12.360704557869983\n",
    "99\t16.294272743826642\n",
    "100\t144.1096440897587''')\n",
    "t1 = pd.read_csv(out, sep='\\\\t')\n",
    "t1.sample(3)"
   ],
   "id": "c839fdb8efc310aa",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(4, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(t1, x='Q', y='Coverage',\n",
    "             color='blue',\n",
    "             ax=ax)\n",
    "ax.set(yscale='log')\n",
    "ax.set_title('H3K27me3 TCell')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "81cdcb2735fd3dcc",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Monocyte H3K27me3 scores",
   "id": "d8e42e8cf1f39f05"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "out = StringIO('''Q\tCoverage\n",
    "1\t0.3674968528850454\n",
    "2\t0.3674968528850454\n",
    "3\t0.3674968528850454\n",
    "4\t0.5783312352566969\n",
    "5\t0.5783312352566969\n",
    "6\t0.5783312352566969\n",
    "7\t0.5783312352566969\n",
    "8\t0.5783312352566969\n",
    "9\t0.5783312352566969\n",
    "10\t0.5783312352566969\n",
    "11\t0.5783312352566969\n",
    "12\t0.5783312352566969\n",
    "13\t0.7891656176283485\n",
    "14\t0.7891656176283485\n",
    "15\t0.7891656176283485\n",
    "16\t0.7891656176283485\n",
    "17\t0.7891656176283485\n",
    "18\t0.7891656176283485\n",
    "19\t0.7891656176283485\n",
    "20\t0.7891656176283485\n",
    "21\t0.7891656176283485\n",
    "22\t0.7891656176283485\n",
    "23\t0.7891656176283485\n",
    "24\t0.7891656176283485\n",
    "25\t0.7891656176283485\n",
    "26\t0.7891656176283485\n",
    "27\t0.7891656176283485\n",
    "28\t0.7891656176283485\n",
    "29\t0.7891656176283485\n",
    "30\t0.7891656176283485\n",
    "31\t0.7891656176283485\n",
    "32\t0.7891656176283485\n",
    "33\t1.0\n",
    "34\t1.0\n",
    "35\t1.0\n",
    "36\t1.0\n",
    "37\t1.0\n",
    "38\t1.0\n",
    "39\t1.0\n",
    "40\t1.0\n",
    "41\t1.0\n",
    "42\t1.0\n",
    "43\t1.0\n",
    "44\t1.0\n",
    "45\t1.0\n",
    "46\t1.0\n",
    "47\t1.0\n",
    "48\t1.0\n",
    "49\t1.0\n",
    "50\t1.0\n",
    "51\t1.0\n",
    "52\t1.0\n",
    "53\t1.0\n",
    "54\t1.0\n",
    "55\t1.0\n",
    "56\t1.0\n",
    "57\t1.0\n",
    "58\t1.0\n",
    "59\t1.0\n",
    "60\t1.0\n",
    "61\t1.0\n",
    "62\t1.0\n",
    "63\t1.3674968528850453\n",
    "64\t1.3674968528850453\n",
    "65\t1.578331235256697\n",
    "66\t1.578331235256697\n",
    "67\t1.578331235256697\n",
    "68\t1.7891656176283486\n",
    "69\t1.7891656176283486\n",
    "70\t1.7891656176283486\n",
    "71\t1.7891656176283486\n",
    "72\t1.7891656176283486\n",
    "73\t1.7891656176283486\n",
    "74\t1.7891656176283486\n",
    "75\t2.0\n",
    "76\t2.0\n",
    "77\t2.0\n",
    "78\t2.0\n",
    "79\t2.0\n",
    "80\t2.0\n",
    "81\t2.0\n",
    "82\t2.0\n",
    "83\t2.0\n",
    "84\t2.156662470513394\n",
    "85\t2.578331235256697\n",
    "86\t2.7891656176283486\n",
    "87\t2.7891656176283486\n",
    "88\t2.7891656176283486\n",
    "89\t3.0\n",
    "90\t3.0\n",
    "91\t3.0\n",
    "92\t3.578331235256697\n",
    "93\t3.7891656176283486\n",
    "94\t4.0\n",
    "95\t4.578331235256697\n",
    "96\t5.0\n",
    "97\t6.367496852885045\n",
    "98\t8.789165617628349\n",
    "99\t14.0\n",
    "100\t68.53735837982705''')\n",
    "t1 = pd.read_csv(out, sep='\\\\t')\n",
    "t1.sample(3)"
   ],
   "id": "72aeb94673028a05",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(4, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(t1, x='Q', y='Coverage',\n",
    "             color='blue',\n",
    "             ax=ax)\n",
    "ax.set(yscale='log')\n",
    "ax.set_title('H3K27me3 Monocytes')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "fa85dd617e8b12fa",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Monocyte CD4 scores",
   "id": "e36b3ac4e8db9ae2"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "out = StringIO('''Q\tCoverage\n",
    "1\t0.8205303096321184\n",
    "2\t1.0171085913600988\n",
    "3\t1.213686873088079\n",
    "4\t1.213686873088079\n",
    "5\t1.213686873088079\n",
    "6\t1.213686873088079\n",
    "7\t1.213686873088079\n",
    "8\t1.4102651548160592\n",
    "9\t1.4102651548160592\n",
    "10\t1.4102651548160592\n",
    "11\t1.4102651548160592\n",
    "12\t1.4102651548160592\n",
    "13\t1.4102651548160592\n",
    "14\t1.4102651548160592\n",
    "15\t1.4102651548160592\n",
    "16\t1.4102651548160592\n",
    "17\t1.6068434365440396\n",
    "18\t1.6068434365440396\n",
    "19\t1.6068434365440396\n",
    "20\t1.6068434365440396\n",
    "21\t1.6068434365440396\n",
    "22\t1.6068434365440396\n",
    "23\t1.6068434365440396\n",
    "24\t1.6068434365440396\n",
    "25\t1.6068434365440396\n",
    "26\t1.6068434365440396\n",
    "27\t1.6068434365440396\n",
    "28\t1.6068434365440396\n",
    "29\t1.6068434365440396\n",
    "30\t1.6068434365440396\n",
    "31\t1.6068434365440396\n",
    "32\t1.8034217182720198\n",
    "33\t1.8034217182720198\n",
    "34\t1.8034217182720198\n",
    "35\t1.8034217182720198\n",
    "36\t1.8034217182720198\n",
    "37\t1.8034217182720198\n",
    "38\t1.8034217182720198\n",
    "39\t1.8034217182720198\n",
    "40\t1.8034217182720198\n",
    "41\t1.8034217182720198\n",
    "42\t1.8034217182720198\n",
    "43\t1.8034217182720198\n",
    "44\t1.8034217182720198\n",
    "45\t1.8034217182720198\n",
    "46\t1.8034217182720198\n",
    "47\t1.8034217182720198\n",
    "48\t1.8034217182720198\n",
    "49\t1.8034217182720198\n",
    "50\t1.8034217182720198\n",
    "51\t2.0\n",
    "52\t2.0\n",
    "53\t2.0\n",
    "54\t2.0\n",
    "55\t2.0\n",
    "56\t2.0\n",
    "57\t2.0\n",
    "58\t2.0\n",
    "59\t2.0\n",
    "60\t2.0\n",
    "61\t2.0\n",
    "62\t2.0\n",
    "63\t2.0\n",
    "64\t2.0\n",
    "65\t2.8205303096321184\n",
    "66\t3.213686873088079\n",
    "67\t3.213686873088079\n",
    "68\t3.4102651548160594\n",
    "69\t3.4102651548160594\n",
    "70\t3.4102651548160594\n",
    "71\t3.4102651548160594\n",
    "72\t3.6068434365440396\n",
    "73\t3.6068434365440396\n",
    "74\t3.6068434365440396\n",
    "75\t3.6068434365440396\n",
    "76\t3.6068434365440396\n",
    "77\t3.6068434365440396\n",
    "78\t3.6068434365440396\n",
    "79\t3.80342171827202\n",
    "80\t3.80342171827202\n",
    "81\t3.80342171827202\n",
    "82\t3.80342171827202\n",
    "83\t3.80342171827202\n",
    "84\t3.80342171827202\n",
    "85\t4.0\n",
    "86\t4.0\n",
    "87\t4.0\n",
    "88\t4.0\n",
    "89\t4.0\n",
    "90\t5.213686873088079\n",
    "91\t5.410265154816059\n",
    "92\t5.60684343654404\n",
    "93\t5.60684343654404\n",
    "94\t5.80342171827202\n",
    "95\t5.80342171827202\n",
    "96\t6.0\n",
    "97\t6.0\n",
    "98\t7.60684343654404\n",
    "99\t8.0\n",
    "100\t245.3337880504351''')\n",
    "t2 = pd.read_csv(out, sep='\\\\t')\n",
    "t2.sample(3)"
   ],
   "id": "a144a4765f0ee304",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(4, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(t2, x='Q', y='Coverage',\n",
    "             color='blue',\n",
    "             ax=ax)\n",
    "ax.set(yscale='log')\n",
    "ax.set_title('H3K27me3 CD4')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "4ea00b28e2c4ad5b",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## From files",
   "id": "808d5fb1dd193e9"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for m in ['CTCF', 'ATAC'] + MODIFICATIONS:\n",
    "    for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "        print(ds, m)\n",
    "        dfdsm = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "                                 (full_sensitivity['modification'] == m) &\n",
    "                                 (full_sensitivity['Gap'] == 0)]\n",
    "        for (cell, rep), t in dfdsm.groupby(['cell', 'replicate']):\n",
    "            if len(t['file'].unique()) != 1:\n",
    "                continue\n",
    "            assert len(t['file'].unique()) == 1\n",
    "            cov = t['file'].values[0].replace('.sensitivity.tsv', '.coverage.tsv')\n",
    "            if not os.path.exists(cov):\n",
    "                continue\n",
    "            t = pd.read_csv(cov, sep='\\t')\n",
    "            t['dataset'] = ds\n",
    "            t['modification'] = m\n",
    "            t['cell'] = cell\n",
    "            t['replicate'] = rep\n",
    "            ts.append(t)\n",
    "cov_df = pd.concat(ts).reset_index(drop=True)\n",
    "del ts\n",
    "cov_df['CoverageI'] = cov_df['Coverage'].astype(int)\n",
    "cov_df.sample(5)"
   ],
   "id": "a04b5749206199d7",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r, ds) in EXAMPLES:\n",
    "    cov = cov_df[(cov_df['dataset'] == ds) &\n",
    "                 (cov_df['modification'] == m) &\n",
    "                 (cov_df['cell'] == c) &\n",
    "                 (cov_df['replicate'] == r)]\n",
    "    plt.figure(figsize=(4, 2))\n",
    "    ax = plt.axes()\n",
    "    sns.lineplot(cov, x='Q', y='CoverageI',\n",
    "                 color='blue',\n",
    "                 ax=ax)\n",
    "    ax.set(yscale='log')\n",
    "    ax.set_title(f'{ds} {m} {c} {r}')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "75f702af4c0ca61c",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## QQs",
   "id": "568dfcf54654c3f1"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "covs = []\n",
    "for (m, c, r, ds) in [\n",
    "    ('H3K27me3', 'Monocyte', '', 'Roadmap'),\n",
    "    ('H3K27me3', 'CD4', '', 'Roadmap'),\n",
    "    ('H3K27me3', 'TCell', '', 'Roadmap')\n",
    "]:\n",
    "    covs.append((f'{ds} {m} {c} {r}',\n",
    "                 cov_df[(cov_df['dataset'] == ds) &\n",
    "                        (cov_df['modification'] == m) &\n",
    "                        (cov_df['cell'] == c) &\n",
    "                        (cov_df['replicate'] == r)]))\n",
    "\n",
    "for i, j in [(0, 1), (0, 2)]:\n",
    "    n1, t1 = covs[i]\n",
    "    n2, t2 = covs[j]\n",
    "    plt.figure(figsize=(3, 3))\n",
    "    ax = plt.axes()\n",
    "    minxy = min(t1['Coverage'].min(), t2['Coverage'].min())\n",
    "    maxxy = max(t1['Coverage'].max(), t2['Coverage'].max())\n",
    "    sns.lineplot(x=[minxy, maxxy], y=[minxy, maxxy], ax=ax, color='black', alpha=0.5, linestyle='dotted')\n",
    "    sns.lineplot(x=t1['Coverage'].tolist(), y=t2['Coverage'].tolist(),\n",
    "                 color='blue',\n",
    "                 sort=False, estimator=None,\n",
    "                 ax=ax)\n",
    "    ax.set(xscale='log')\n",
    "    ax.set(yscale='log')\n",
    "    ax.set_xlim(minxy, maxxy)\n",
    "    ax.set_ylim(minxy, maxxy)\n",
    "    ax.set_xlabel(n1)\n",
    "    ax.set_ylabel(n2)\n",
    "    ax.set_title('QQ')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "d449c02ec74786b0",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## All",
   "id": "d0d52a2b5a60ddb1"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(5, 2))\n",
    "ax = plt.axes()\n",
    "sns.lineplot(cov_df[cov_df['Q'] >= 50],\n",
    "             x='Q', y='CoverageI',\n",
    "             hue='modification',\n",
    "             hue_order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "             alpha=0.5, errorbar=('se', .95),\n",
    "             ax=ax)\n",
    "ax.set_title('Coverage')\n",
    "ax.set(yscale='log')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "id": "eb5736e7f20b78f0",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.rc('font', size=5)\n",
    "plt.figure(figsize=(3, 6))\n",
    "axs = [plt.subplot(7, 1, i + 1) for i in range(7)]\n",
    "for i, m in enumerate(['CTCF', 'ATAC'] + MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    sns.lineplot(cov_df[cov_df['modification'] == m],\n",
    "                 x='Q', y='CoverageI',\n",
    "                 hue='dataset',\n",
    "                 alpha=0.5, errorbar=('se', .95),\n",
    "                 ax=ax)\n",
    "    ax.set_title(f'{m} Coverage')\n",
    "    ax.set(yscale='log')\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "plt.rc('font', size=8)"
   ],
   "id": "2d96aa7df14f7a86",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Candidates from bw coverage",
   "id": "5e5dfb90c017211"
  },
  {
   "metadata": {},
   "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"
   ],
   "id": "dc05d112fa85ed1b",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "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)"
   ],
   "id": "32c7062ac6dc7c2c",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "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)"
   ],
   "id": "381b59d795023044",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Prepare BAM files\n",
    "# https://github.com/Boyle-Lab/Blacklist/raw/refs/heads/master/lists/hg38-blacklist.v2.bed.gz\n",
    "BLACKLIST_PATH = os.path.expanduser('~/data') + '/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"
   ],
   "id": "abf8741868196421",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def collect_bw_stats(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(\n",
    "                    [x or 0 for x in bw.stats(chr, 1, CHROM_SIZES[chr], nBins=iend - istart, exact=True, type='sum')])\n",
    "\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",
    "    region_stats = []\n",
    "    with pyBigWig.open(bw_path) as bw:\n",
    "        for _, (chr, start, end) in df.iterrows():\n",
    "            stats = np.array(\n",
    "                bw.stats(chr, start, end, nBins=int(math.ceil((end - start) / resolution)), exact=True, type='sum'))\n",
    "            region_stats.append(stats)\n",
    "    return df, region_stats"
   ],
   "id": "58f529620ddb8873",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "REGION_LEN = 20_000\n",
    "TOP_REGIONS = 10_000\n",
    "WORK_REGIONS = 1000\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)}"
   ],
   "id": "ab73a293aec02592",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Candidates calling from ModelToPeaks.kt",
   "id": "27f34a2fb4ac442d"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def next_clear_bit(arr, start):\n",
    "    for i in range(start, len(arr)):\n",
    "        if not arr[i]:\n",
    "            return i\n",
    "    return -1\n",
    "\n",
    "\n",
    "def next_set_bit(arr, start):\n",
    "    for i in range(start, len(arr)):\n",
    "        if arr[i]:\n",
    "            return i\n",
    "    return -1\n",
    "\n",
    "\n",
    "def set_bits(arr, start, end):\n",
    "    for i in range(start, end):\n",
    "        arr[i] = True\n",
    "\n",
    "\n",
    "# Copy of BitList.aggregate\n",
    "def aggregate(arr, gap=0):\n",
    "    ranges = []\n",
    "    offset = 0\n",
    "    while offset < len(arr):\n",
    "        left = next_set_bit(arr, offset)\n",
    "        if left == -1:\n",
    "            break\n",
    "        right = min(next_clear_bit(arr, left + 1), len(arr))\n",
    "        if right == -1:\n",
    "            break\n",
    "        while gap > 0 and right < len(arr):\n",
    "            next_set = next_set_bit(arr, right + 1)\n",
    "            if next_set != -1 and next_set - right <= gap:\n",
    "                right = min(next_clear_bit(arr, next_set + 1), len(arr))\n",
    "            else:\n",
    "                break\n",
    "        ranges.append((left, right))\n",
    "        offset = right\n",
    "    return ranges"
   ],
   "id": "a382dec873662ae8",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def call_candidates(bins, gap):\n",
    "    candidates = aggregate(bins)\n",
    "    if len(candidates) == 0:\n",
    "        return candidates\n",
    "    last_end = -1\n",
    "    for (start, end) in candidates:\n",
    "        if last_end != -1:\n",
    "            if start - last_end < gap:\n",
    "                set_bits(bins, last_end, start)\n",
    "        last_end = end\n",
    "\n",
    "    candidates = aggregate(bins)\n",
    "    for (start, end) in candidates:\n",
    "        assert start < end, f'{start} >= {end}\\n{bins}'\n",
    "    return candidates\n",
    "\n",
    "\n",
    "def analyze_candidates(coverages, fdr=0.05):\n",
    "    rows = []\n",
    "    bins = np.ndarray(len(coverages), dtype='bool')\n",
    "    for s in np.exp(np.linspace(-8, 1, 100) * math.log(10)):\n",
    "        cq = np.percentile(coverages, 100 * (1.0 - s * fdr))\n",
    "        for gap in [0]:\n",
    "            for i in range(len(bins)):\n",
    "                bins[i] = coverages[i] >= cq\n",
    "            candidates = call_candidates(bins, gap)\n",
    "            avg_len = sum((end - start) for (start, end) in candidates) / len(candidates) if len(candidates) > 0 else 0\n",
    "            rows.append((s, gap, len(candidates), avg_len))\n",
    "    df = pd.DataFrame(rows, columns=['Sensitivity', 'Gap', 'CandidatesN', 'CandidatesAL'])\n",
    "    return df"
   ],
   "id": "b18db0b62e28eb8d",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Visualization",
   "id": "6465e47fd08c4cb4"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def plot_candidates(t, title=None):\n",
    "    plt.figure(figsize=(5, 4))\n",
    "    ax = plt.axes()\n",
    "    sns.lineplot(data=t, x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap', estimator=None,\n",
    "                 # marker='o', markersize=3,\n",
    "                 palette='tab20',\n",
    "                 sort=False,\n",
    "                 ax=ax)\n",
    "    ax.set(yscale='log')\n",
    "    ax.set(xscale='log')\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    if title:\n",
    "        plt.title(title)\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "2fad75334b4dfc0f",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Examples",
   "id": "d3e9879e6ff01306"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r, ds) in EXAMPLES:\n",
    "    bw_path = all_bws[(all_bws['modification'] == m) &\n",
    "                      (all_bws['cell'] == c) &\n",
    "                      (all_bws['replicate'] == r) &\n",
    "                      (all_bws['dataset'] == ds)]['file'].values[0]\n",
    "\n",
    "    _, region_stats = collect_bw_stats(bw_path, chrom_sizes_part, blacklist_df, REGION_LEN, TOP_REGIONS,\n",
    "                                       WORK_REGIONS, RESOLUTION)\n",
    "    coverages = np.concatenate(region_stats, axis=0)\n",
    "    t = analyze_candidates(coverages)\n",
    "    display(t.head(10))\n",
    "    display(t.tail(20))\n",
    "    # For log\n",
    "    t.fillna(0, inplace=True)\n",
    "    t['CandidatesN'] += 1\n",
    "    t['CandidatesAL'] += 1\n",
    "    # break\n",
    "    t.sort_values(['Sensitivity', 'Gap'], inplace=True)\n",
    "    plot_candidates(t, f'{m} {c} {r}')\n",
    "    plot_projections(t, f'{m} {c} {r}', None)\n"
   ],
   "id": "e3b369875e793994",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Segments",
   "id": "99f1e20f2ce61ce5"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r, ds) in EXAMPLES:\n",
    "    bw_path = all_bws[(all_bws['modification'] == m) &\n",
    "                      (all_bws['cell'] == c) &\n",
    "                      (all_bws['replicate'] == r) &\n",
    "                      (all_bws['dataset'] == ds)]['file'].values[0]\n",
    "\n",
    "    _, region_stats = collect_bw_stats(bw_path, chrom_sizes_part, blacklist_df, REGION_LEN, TOP_REGIONS,\n",
    "                                       WORK_REGIONS, RESOLUTION)\n",
    "    coverages = np.concatenate(region_stats, axis=0)\n",
    "    rows = []\n",
    "    fdr = 0.05\n",
    "    prev_bins = np.ndarray(len(coverages), dtype='bool')\n",
    "    bins = np.ndarray(len(coverages), dtype='bool')\n",
    "    for i, s in enumerate(np.exp(np.linspace(-8, 1, 50) * math.log(10))):\n",
    "        cq = np.percentile(coverages, 100 * (1.0 - s * fdr))\n",
    "        if i == 0:\n",
    "            for j in range(len(bins)):\n",
    "                prev_bins[j] = coverages[j] >= cq\n",
    "            continue\n",
    "        for j in range(len(bins)):\n",
    "            bins[j] = coverages[j] >= cq\n",
    "        if prev_bins is None:\n",
    "            prev_bins = bins.copy()\n",
    "            continue\n",
    "        candidates = call_candidates(bins, 0)\n",
    "        total = len(candidates)\n",
    "        overlap_candidates = 0\n",
    "        for (start, end) in candidates:\n",
    "            ovlp = any(prev_bins[i] for i in range(start, end))\n",
    "            if ovlp:\n",
    "                overlap_candidates += 1\n",
    "        rows.append((s, total, total - overlap_candidates, overlap_candidates))\n",
    "        for j in range(len(bins)):\n",
    "            prev_bins[j] = bins[j]\n",
    "\n",
    "    sdf = pd.DataFrame(rows, columns=['Sensitivity', 'Total', 'New', 'Old'])\n",
    "    # display(sdf)\n",
    "    sdf['New %'] = [d(n, t) * 100 for n, t in zip(sdf['New'], sdf['Total'])]\n",
    "    sdf['Total %'] = 100\n",
    "    plt.figure(figsize=(10, 1))\n",
    "    ax = plt.axes()\n",
    "    # Background\n",
    "    sns.barplot(sdf, x='Sensitivity', y='Total %',\n",
    "                color='orange',\n",
    "                ax=ax)\n",
    "    # Foreground\n",
    "    sns.barplot(sdf, x='Sensitivity', y='New %',\n",
    "                color='blue',\n",
    "                ax=ax)\n",
    "    ax.set_ylabel('peaks %')\n",
    "    name = f'{ds} {m} {c} {r}'\n",
    "    ax.set_title(name)\n",
    "    # Points\n",
    "    # sp = sps[name]\n",
    "    # ax.axvline(x=sp[2], ymin=0, ymax=100, color='blue', lw=2)\n",
    "    # ax.axvline(x=sp[1], ymin=0, ymax=100, color='red', lw=2)\n",
    "    # ax.axvline(x=sp[0], ymin=0, ymax=100, color='green', lw=2)\n",
    "    # Limit ticks number\n",
    "    for i, label in enumerate(ax.get_xticklabels()):\n",
    "        if 0 < i < len(sdf) - 1 and i % 20 != 0:\n",
    "            label.set_visible(False)\n",
    "    ax.invert_xaxis()\n",
    "    # add legend\n",
    "    top_bar = mpatches.Patch(color='orange', label='old %')\n",
    "    bottom_bar = mpatches.Patch(color='blue', label='new %')\n",
    "    plt.legend(handles=[bottom_bar, top_bar])\n",
    "    sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.show()\n"
   ],
   "id": "9c362a077858542e",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "2818d2c7f44c2431",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "2cc828d0dbd7ee9f",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "##  All",
   "id": "c65f4f30fafb8ae"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# for m in ['ATAC'] + MODIFICATIONS:\n",
    "#     for ds in ['ENCODE', 'Roadmap', 'ImmGen']:\n",
    "#         print(ds, m)\n",
    "#         dft = all_bws[(all_bws['dataset'] == ds) & (all_bws['modification'] == m)].copy().reset_index(drop=True)\n",
    "#\n",
    "#         n = len(dft)\n",
    "#         print(n)\n",
    "#         if n <= 6:\n",
    "#             ncols, nrows = n, 1\n",
    "#         else:\n",
    "#             ncols = int(math.floor(math.sqrt(n) + 4))\n",
    "#             nrows = int(math.ceil(n / ncols))\n",
    "#         plt.figure(figsize=(3 * ncols + 1, 3 * nrows))\n",
    "#         axs = [plt.subplot(nrows, ncols, i + 1) for i in range(n)]\n",
    "#         for i, dftr in dft.iterrows():\n",
    "#             c, r, bw_path = dftr['cell'], dftr['replicate'], dftr['file']\n",
    "#             print(ds, m, c, r, bw_path)\n",
    "#             regions_df, region_stats = collect_bw_stats(\n",
    "#                 bw_path, chrom_sizes_part, blacklist_df, REGION_LEN, TOP_REGIONS, WORK_REGIONS, RESOLUTION\n",
    "#             )\n",
    "#             coverages = np.concatenate(region_stats, axis=0)\n",
    "#             tt = analyze_candidates(coverages)\n",
    "#             ax = axs[i]\n",
    "#             sns.lineplot(data=tt, x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap', estimator=None,\n",
    "#                      marker='o', markersize=3,\n",
    "#                      palette='tab20',\n",
    "#                     sort=False,\n",
    "#                     ax=ax)\n",
    "#             for x, y, s in zip(tt['CandidatesN'], tt['CandidatesAL'], tt['Sensitivity']):\n",
    "#                 if s in [1e-3, 0.1]:\n",
    "#                     ax.text(x, y, str(s), fontsize=5)\n",
    "#             # pivot = full_sensitivity[(full_sensitivity['dataset'] == ds) &\n",
    "#             #                   (full_sensitivity['modification'] == m) &\n",
    "#             #                   (full_sensitivity['cell'] == c) &\n",
    "#             #                   (full_sensitivity['replicate'] == r)]['pivot'].values[0]\n",
    "#             # ax.set_title(f'{ds} {m} {c} {r} {\"\" if pivot == \"null\" else \"\"}')\n",
    "#             ax.set_title(f'{ds} {m} {c} {r}')\n",
    "#             ax.set(yscale='log')\n",
    "#             ax.set(xscale='log')\n",
    "#             if i == n - 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()"
   ],
   "id": "531af916bb9e0ad1",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Chips quality analysis",
   "id": "b092045df7cc4341"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "CHIPS_PATH = os.path.expanduser('~/data/2023_chips')\n",
    "N = 5\n",
    "MULTIPLIERS = [1.0, 0.7, 0.5, 0.2, 0.1]"
   ],
   "id": "9ec64cdf49eda85b",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def load_chips_sensitivity(path):\n",
    "    ts = []\n",
    "    for (m, n, mult) in product(\n",
    "            ['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3', 'mixed'],\n",
    "            range(1, N + 1),\n",
    "            MULTIPLIERS):\n",
    "        res = glob.glob(f'{path}/{m}*_{n}_{mult}_*.sensitivity.tsv')\n",
    "        if len(res) == 0:\n",
    "            print(f'Nothing found for {path}/{m}*_{n}_{mult}_*.sensitivity.tsv')\n",
    "            continue\n",
    "        if len(res) > 1:\n",
    "            print(f'More than 1 found for {path}/{m}*_{n}_{mult}_*.sensitivity.tsv:\\n{res}')\n",
    "            continue\n",
    "        tsv_path = res[0]\n",
    "        print(tsv_path, m, n, mult)\n",
    "        t = pd.read_csv(tsv_path, sep='\\t')\n",
    "        t['file'] = tsv_path\n",
    "        t['modification'] = m\n",
    "        t['n'] = n\n",
    "        t['multiplier'] = mult\n",
    "        ts.append(t)\n",
    "    df = pd.concat(ts).reset_index(drop=True)\n",
    "    return df"
   ],
   "id": "b558dd11bc4f1f99",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "chips_sens = load_chips_sensitivity(CHIPS_PATH + '/span')\n",
    "# For log\n",
    "chips_sens['CandidatesN'] += 1\n",
    "chips_sens['CandidatesAL'] += 1\n",
    "chips_sens.sample(3)"
   ],
   "id": "dfd822eca40b2960",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for m in MODIFICATIONS:\n",
    "    print(m)\n",
    "    dfdsm = chips_sens[chips_sens['modification'] == m].copy()\n",
    "    dfdsm.sort_values(by=['multiplier', 'Sensitivity'], ascending=False, inplace=True)\n",
    "    # t = dfdsm[(dfdsm['Sensitivity'] >= 1e-3) & (dfdsm['Sensitivity']<=1)]\n",
    "    t = dfdsm\n",
    "    minx, maxx = t['CandidatesN'].min(), t['CandidatesN'].max()\n",
    "    miny, maxy = t['CandidatesAL'].min(), t['CandidatesAL'].max()\n",
    "    n = len(dfdsm[['n', 'multiplier']].copy().drop_duplicates())\n",
    "    print(n)\n",
    "    if n <= 6:\n",
    "        ncols, nrows = n, 1\n",
    "    else:\n",
    "        ncols = int(math.floor(math.sqrt(n) + 4))\n",
    "        nrows = int(math.ceil(n / ncols))\n",
    "    plt.figure(figsize=(3 * ncols + 1, 3 * nrows))\n",
    "    axs = [plt.subplot(nrows, ncols, i + 1) for i in range(n)]\n",
    "    for i, ((mult, n), t) in enumerate(dfdsm.groupby(['multiplier', 'n'])):\n",
    "        # tt = t[(t['Sensitivity'] >= 1e-3) & (t['Sensitivity']<=1)].copy()\n",
    "        tt = t\n",
    "        ax = axs[i]\n",
    "        sns.lineplot(data=tt, x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap',\n",
    "                     # marker='o', markersize=2,\n",
    "                     palette='tab20',\n",
    "                     sort=False,\n",
    "                     alpha=0.5,\n",
    "                     ax=ax)\n",
    "        # sns.scatterplot(data=tt,\n",
    "        #                 x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap', s=10, palette='tab20',\n",
    "        #                 legend=False, ax=ax)\n",
    "        # sns.scatterplot(data=tt,\n",
    "        #                 x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap', s=30, palette='tab20',\n",
    "        #                 legend=False, ax=ax)\n",
    "        # for x, y, v in zip(tt['CandidatesN'], tt['CandidatesAL'], tt['Sensitivity']):\n",
    "        #     if v in [1.0, 0.001]:\n",
    "        #         ax.text(x, y, str(v), fontsize=5)\n",
    "        ax.set(yscale='log')\n",
    "        ax.set(xscale='log')\n",
    "        # pivot = tt['pivot'].unique()[0]\n",
    "        # ax.set_title(f'{m} {mult} {n} {f\"P {pivot}\" if pivot != \"null\" else \"\"}')\n",
    "        ax.set_title(f'{m} {mult} {n}')\n",
    "        ax.set_xlim(minx, maxx)\n",
    "        ax.set_ylim(miny, maxy)\n",
    "        if i == n - 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()"
   ],
   "id": "8bd552875946d4ad",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, n) in [\n",
    "    ('H3K4me3', 1),\n",
    "]:\n",
    "    t = chips_sens[(chips_sens['modification'] == m) &\n",
    "                   (chips_sens['n'] == n)].copy()\n",
    "    t.sort_values(by=['Sensitivity'], ascending=False, inplace=True)\n",
    "    plt.figure(figsize=(15, 3))\n",
    "    axs = [plt.subplot(1, len(MULTIPLIERS), i + 1) for i in range(len(MULTIPLIERS))]\n",
    "    for i, mult in enumerate(sorted(MULTIPLIERS, reverse=True)):\n",
    "        ax = axs[i]\n",
    "        tt = t[t['multiplier'] == mult]\n",
    "        # tt = t[(t['Sensitivity'] >= 1e-3) & (t['Sensitivity'] <= 1.0)]\n",
    "        # tt = t\n",
    "        sns.lineplot(data=tt, x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap', estimator=None,\n",
    "                     # marker='o', markersize=3,\n",
    "                     palette='tab20',\n",
    "                     sort=False,\n",
    "                     alpha=0.5,\n",
    "                     ax=ax)\n",
    "        ttt = tt[tt['Gap'] == 0].sort_values(by=['Sensitivity']).copy().reset_index(drop=True)\n",
    "        sp = detect_sensitivity_triangle(ttt)\n",
    "        sns.scatterplot(data=ttt[ttt.index == sp[0]],\n",
    "                        x=\"CandidatesN\", y=\"CandidatesAL\", color='green', s=20,\n",
    "                        legend=False, ax=ax)\n",
    "        sns.scatterplot(data=ttt[ttt.index == sp[1]],\n",
    "                        x=\"CandidatesN\", y=\"CandidatesAL\", color='red', s=40,\n",
    "                        legend=False, ax=ax)\n",
    "        sns.scatterplot(data=ttt[ttt.index == sp[2]],\n",
    "                        x=\"CandidatesN\", y=\"CandidatesAL\", color='blue', s=20,\n",
    "                        legend=False, ax=ax)\n",
    "\n",
    "        # sns.scatterplot(data=tt[tt['Sensitivity'].isin([1.0, 0.001])],\n",
    "        #                 x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap', s=10, palette='tab20',\n",
    "        #                 legend=False, ax=ax)\n",
    "        # sns.scatterplot(data=tt[tt['Sensitivity'].isin([0.1])],\n",
    "        #                 x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap', s=30, palette='tab20',\n",
    "        #                 legend=False, ax=ax)\n",
    "        # for x, y, v in zip(tt['CandidatesN'], tt['CandidatesAL'], tt['Sensitivity']):\n",
    "        #     if v in [1.0, 0.001]:\n",
    "        #         ax.text(x, y, str(v), fontsize=5)\n",
    "        ax.set(yscale='log')\n",
    "        ax.set(xscale='log')\n",
    "        # pivot = tt['pivot'].values[0]\n",
    "        # ax.set_title(f'{m} {mult} {n} {f\"P {pivot}\" if pivot != \"null\" else \"\"}')\n",
    "        ax.set_title(f'{m} {mult} {n}')\n",
    "        if i == len(MULTIPLIERS) - 1:\n",
    "            sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "        else:\n",
    "            ax.legend().set_visible(False)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "fd28b6babfe2041e",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Mix control",
   "id": "236f775db0a18fd5"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "MIX_PATH = os.path.expanduser('~/data/2023_mix_chr1')\n",
    "MIX_SIGNAL_FRACTIONS = [0.7, 0.5, 0.2, 0.1]"
   ],
   "id": "2b6ef215490c6128",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def load_mix_sensitivity(path):\n",
    "    ts = []\n",
    "    for (m, c, r, f) in product(\n",
    "            ['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3'],\n",
    "            IMMUNE_CELLS,\n",
    "            IMMUNE_REPS,\n",
    "            MIX_SIGNAL_FRACTIONS,\n",
    "    ):\n",
    "        res = glob.glob(f'{path}/{m}_{c}_{r}*_{f}_*.sensitivity.tsv')\n",
    "        if r == '':\n",
    "            res = [f for f in res if len([ro for ro in IMMUNE_REPS if ro != r and ro in f]) == 0]\n",
    "        if len(res) == 0:\n",
    "            print(f'Nothing found for {path}/{m}_{c}_{r}*_{f}_*.sensitivity.tsv')\n",
    "            continue\n",
    "        if len(res) > 1:\n",
    "            print(f'More than 1 found for {path}/{m}_{c}_{r}*_{f}_*.sensitivity.tsv:\\n{res}')\n",
    "            continue\n",
    "        tsv_path = res[0]\n",
    "        print(tsv_path, m, f)\n",
    "        t = pd.read_csv(tsv_path, sep='\\t')\n",
    "        t['file'] = tsv_path\n",
    "        t['modification'] = m\n",
    "        t['cell'] = c\n",
    "        t['replicate'] = r\n",
    "        t['fraction'] = f\n",
    "        t['SensitivityN'] = range(len(t))\n",
    "        ts.append(t)\n",
    "    df = pd.concat(ts).reset_index(drop=True)\n",
    "    return df"
   ],
   "id": "fad22c6f9731f8b4",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "mix_sens = load_mix_sensitivity(MIX_PATH + '/span')\n",
    "# For log\n",
    "mix_sens.fillna(0, inplace=True)\n",
    "mix_sens['CandidatesN'] += 1\n",
    "mix_sens['CandidatesAL'] += 1\n",
    "mix_sens.sample(3)"
   ],
   "id": "ec2f61dc011501b",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Examples\n",
    "sps_mix = {}\n",
    "for (m, c, r) in [\n",
    "    ('H3K4me3', 'BCell', ''),\n",
    "    # ('H3K4me3', 'CD4', ''),\n",
    "    ('H3K27ac', 'CD34', ''),\n",
    "    ('H3K4me1', 'CD4', 'rep1'),\n",
    "    ('H3K27me3', 'NK', ''),\n",
    "    ('H3K36me3', 'BCell', ''),\n",
    "]:\n",
    "    t = mix_sens[(mix_sens['modification'] == m) &\n",
    "                 (mix_sens['cell'] == c) &\n",
    "                 (mix_sens['replicate'] == r)].copy()\n",
    "    t.sort_values(by=['fraction', 'Sensitivity'], ascending=False, inplace=True)\n",
    "    plt.figure(figsize=(12, 3))\n",
    "    axs = [plt.subplot(1, 4, i + 1) for i in range(4)]\n",
    "    for i, f in enumerate(sorted(MIX_SIGNAL_FRACTIONS, reverse=True)):\n",
    "        ax = axs[i]\n",
    "        tt = t[t['fraction'] == f]\n",
    "        assert len(tt['file'].unique()) == 1\n",
    "        sns.lineplot(data=tt, x=\"CandidatesN\", y=\"CandidatesAL\", hue='Gap', estimator=None,\n",
    "                     palette='tab20',\n",
    "                     sort=False,\n",
    "                     alpha=0.5,\n",
    "                     ax=ax)\n",
    "        ttt = tt[tt['Gap'] == 0].sort_values(by=['Sensitivity']).copy().reset_index(drop=True)\n",
    "        sp = detect_sensitivity_triangle(ttt)\n",
    "        sps_mix[(m, c, r, i)] = (tt['file'].values[0], ttt, sp)\n",
    "        sns.scatterplot(data=ttt[ttt.index == sp[0]],\n",
    "                        x=\"CandidatesN\", y=\"CandidatesAL\", color='green', s=20,\n",
    "                        legend=False, ax=ax)\n",
    "        sns.scatterplot(data=ttt[ttt.index == sp[1]],\n",
    "                        x=\"CandidatesN\", y=\"CandidatesAL\", color='red', s=40,\n",
    "                        legend=False, ax=ax)\n",
    "        sns.scatterplot(data=ttt[ttt.index == sp[2]],\n",
    "                        x=\"CandidatesN\", y=\"CandidatesAL\", color='blue', s=20,\n",
    "                        legend=False, ax=ax)\n",
    "        ax.set(yscale='log')\n",
    "        ax.set(xscale='log')\n",
    "        # pivot = tt['pivot'].unique()[0]\n",
    "        # ax.set_title(f'{m} {c} {r} {f} {f\"P {pivot}\" if pivot != \"null\" else \"\"}')\n",
    "        ax.set_title(f'{m} {c} {r} {f}')\n",
    "        if i == len(axs) - 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()"
   ],
   "id": "9cea991e5e7626c7",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "mfs = list(sorted(MIX_SIGNAL_FRACTIONS, reverse=True))\n",
    "ts = []\n",
    "for (m, c, r, i), (file, df, sp) in sps_mix.items():\n",
    "    print(m, c, r, mfs[i])\n",
    "    df['modification'] = m\n",
    "    df['cell'] = c\n",
    "    df['replicate'] = r\n",
    "    ts.append(df)\n",
    "sps_mix_df = pd.concat(ts).reset_index(drop=True)\n",
    "del ts"
   ],
   "id": "86f834db9ed303e7",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Examples\n",
    "for (m, c, r) in [\n",
    "    ('H3K4me3', 'BCell', ''),\n",
    "    # ('H3K4me3', 'CD4', ''),\n",
    "    ('H3K27ac', 'CD34', ''),\n",
    "    ('H3K4me1', 'CD4', 'rep1'),\n",
    "    ('H3K27me3', 'NK', ''),\n",
    "    ('H3K36me3', 'BCell', ''),\n",
    "]:\n",
    "    plt.figure(figsize=(15, 1.5))\n",
    "    axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "    t = sps_mix_df[(sps_mix_df['modification'] == m) & (sps_mix_df['cell'] == c) & (sps_mix_df['replicate'] == r)]\n",
    "    for i, y in enumerate(['CandidatesN', 'CandidatesAL', 'CandidatesML', 'SignalNoiseRatio', 'SignalControlRatio']):\n",
    "        ax = axs[i]\n",
    "        sns.lineplot(data=t, x='SensitivityN', y=y, estimator=None,\n",
    "                     hue='fraction',\n",
    "                     hue_order=[0.7, 0.5, 0.2, 0.1],\n",
    "                     # palette='tab20',\n",
    "                     alpha=0.8,\n",
    "                     sort=False,\n",
    "                     ax=ax)\n",
    "        # ax.set(xscale='log')\n",
    "        ax.set(yscale='log')\n",
    "        # ax.invert_xaxis()\n",
    "        ax.set_title(f'{m} {c} {r}')\n",
    "        if i == len(axs) - 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()"
   ],
   "id": "7873e6a965b3b6b6",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for (m, c, r) in [\n",
    "    ('H3K4me3', 'BCell', ''),\n",
    "    ('H3K27ac', 'CD34', ''),\n",
    "    ('H3K4me1', 'CD4', 'rep1'),\n",
    "    ('H3K27me3', 'NK', ''),\n",
    "    ('H3K36me3', 'BCell', ''),\n",
    "]:\n",
    "    plt.figure(figsize=(20, 2))\n",
    "    axs = [plt.subplot(1, 4, i + 1) for i in range(4)]\n",
    "    for i, f in enumerate(sorted(MIX_SIGNAL_FRACTIONS, reverse=True)):\n",
    "        ax = axs[i]\n",
    "        file, _, sp = sps_mix[(m, c, r, i)]\n",
    "        segments = file.replace('.sensitivity.tsv', '.segments.tsv')\n",
    "        sdf = pd.read_csv(segments, sep='\\t')\n",
    "        sdf['New %'] = [d(n, t) * 100 for n, t in zip(sdf['New'], sdf['CandidatesN'])]\n",
    "        sdf['Total %'] = 100\n",
    "        sdf['SensitivityN'] = range(len(sdf))\n",
    "        # Background\n",
    "        sns.barplot(sdf, x='SensitivityN', y='Total %',\n",
    "                    color='orange',\n",
    "                    ax=ax)\n",
    "        # Foreground\n",
    "        sns.barplot(sdf, x='SensitivityN', y='New %',\n",
    "                    color='blue',\n",
    "                    ax=ax)\n",
    "        ax.set_ylabel('peaks %')\n",
    "        name = f'{ds} {m} {c} {r} {f}'\n",
    "        ax.set_title(name)\n",
    "        # Points\n",
    "        ax.axvline(x=sp[2], ymin=0, ymax=100, color='blue', lw=2)\n",
    "        ax.axvline(x=sp[1], ymin=0, ymax=100, color='red', lw=2)\n",
    "        ax.axvline(x=sp[0], ymin=0, ymax=100, color='green', lw=2)\n",
    "        # Limit ticks number\n",
    "        for i, label in enumerate(ax.get_xticklabels()):\n",
    "            label.set_visible(False)\n",
    "        # add legend\n",
    "        if i == len(axs) - 1:\n",
    "            top_bar = mpatches.Patch(color='orange', label='old %')\n",
    "            bottom_bar = mpatches.Patch(color='blue', label='new %')\n",
    "            plt.legend(handles=[bottom_bar, top_bar])\n",
    "            sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "        new_value = 1\n",
    "        for patch in ax.patches:\n",
    "            current_width = patch.get_width()\n",
    "            diff = current_width - new_value\n",
    "            # we change the bar width\n",
    "            patch.set_width(new_value)\n",
    "            # we recenter the bar\n",
    "            patch.set_x(patch.get_x() + diff * .5)\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "8d7c59e6d83017fe",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "3c0b4de7e13d0aa0",
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}