2023 - noise_mix.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noise mix\n",
    "\n",
    "See `~/data/2023_mix_chr1` folder."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'\n",
    "\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from IPython.display import display\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "sns.set_style(\"whitegrid\")\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import re\n",
    "from itertools import product\n",
    "import tempfile"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "def bedl(file):\n",
    "    try:\n",
    "        tf = pd.read_csv(file, sep='\\t', header=None)\n",
    "        return tf[2] - tf[1]\n",
    "    except:\n",
    "        return np.zeros(0)  # Empty file\n",
    "\n",
    "\n",
    "def lines(file):\n",
    "    try:\n",
    "        tf = pd.read_csv(file, sep='\\t', header=None)\n",
    "        return len(tf)\n",
    "    except:\n",
    "        return 0  # Empty file\n",
    "\n",
    "\n",
    "def d(a, b):\n",
    "    return a / b if b != 0 else 0\n",
    "\n",
    "\n",
    "def last_col(file):\n",
    "    try:\n",
    "        cols = len(pd.read_csv(file, sep='\\t', nrows=1, header=None).columns)\n",
    "        return pd.read_csv(file, sep='\\t', header=None, usecols=[cols - 1])[cols - 1]\n",
    "    except:\n",
    "        return np.zeros(0)  # Empty file\n",
    "\n",
    "\n",
    "def sorted_file(file):\n",
    "    ts = tempfile.mktemp()\n",
    "    ! cat {file} | sort -k1,1 -k2,2n -k3,3n > {ts}\n",
    "    return ts"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Load hg38 peaks\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "GSE26320_PATH = os.path.expanduser('~/data/2023_GSE26320')\n",
    "GSE26320_CELLS = ['GM12878', 'HMEC', 'HSMM', 'K562', 'NHEK', 'NHLF', 'H1', 'Huvec', 'HepG2']\n",
    "MODIFICATIONS = ['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']\n",
    "GSE26320_REPS = ['rep1', 'rep2']\n",
    "\n",
    "IMMUNE_PATH = os.path.expanduser('~/data/2023_Immune')\n",
    "IMMUNE_CELLS = ['CD4ABT', 'TCellBB', 'BCell', 'TCell', 'Monocyte', 'PBMC', 'NK', 'CD34', 'CD4', ]  ## Longest first\n",
    "IMMUNE_REPS = ['rep1', 'rep2', 'rep3', '']\n",
    "\n",
    "! mkdir -p {GSE26320_PATH}/pics\n",
    "! mkdir -p {IMMUNE_PATH}/pics"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "TOOLS = ['MACS2', 'MACS2 broad', 'SPAN', 'SICER']\n",
    "palette = plt.cm.get_cmap('tab10')\n",
    "TOOLS_PALETTE = {t: palette(i) for i, t in enumerate(TOOLS)}"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "def load_peaks(path, suffix, modifications, cells, replicates):\n",
    "    df = pd.DataFrame(columns=['file', 'modification', 'cell', 'replicate', 'peaks'],\n",
    "                      dtype=object)\n",
    "    for f in tqdm(os.listdir(path)):\n",
    "        if not f.endswith(suffix):\n",
    "            continue\n",
    "        cell = next((c for c in cells if c in f), None)\n",
    "        mod = next((m for m in modifications if m in f), None)\n",
    "        rep = next((r for r in replicates if r in f), None)\n",
    "        if cell and mod and rep is not None:\n",
    "            peaks_path = os.path.join(path, f)\n",
    "            ps = lines(peaks_path)\n",
    "            df.loc[len(df)] = (peaks_path, mod, cell, rep, ps)\n",
    "    return df"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "def load_peaks_path(path, modifications, cells, replicates):\n",
    "    df_macs2 = load_peaks(os.path.join(path, 'macs2'), '.narrowPeak', modifications, cells, replicates)\n",
    "    df_macs2['tool'] = 'MACS2'\n",
    "    print('MACS2', len(df_macs2))\n",
    "\n",
    "    df_macs2broad = load_peaks(os.path.join(path, 'macs2'), '.broadPeak', modifications, cells, replicates)\n",
    "    df_macs2broad['tool'] = 'MACS2 broad'\n",
    "    print('MACS2 broad', len(df_macs2broad))\n",
    "\n",
    "    df_sicer = load_peaks(os.path.join(path, 'sicer'), 'summary-FDR0.01', modifications, cells, replicates)\n",
    "    df_sicer['tool'] = 'SICER'\n",
    "    print('SICER', len(df_sicer))\n",
    "\n",
    "    df_span = load_peaks(os.path.join(path, 'span'), '.peak', modifications, cells, replicates)\n",
    "    df_span['tool'] = 'SPAN'\n",
    "    print('SPAN', len(df_span))\n",
    "\n",
    "    return pd.concat([df_macs2, df_macs2broad, df_sicer, df_span]).reset_index(drop=True)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_encode_peaks = load_peaks_path(GSE26320_PATH, MODIFICATIONS, GSE26320_CELLS, GSE26320_REPS)\n",
    "df_encode_peaks['dataset'] = 'ENCODE'\n",
    "df_encode_peaks.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_immune_peaks = load_peaks_path(IMMUNE_PATH, MODIFICATIONS, IMMUNE_CELLS, IMMUNE_REPS)\n",
    "df_immune_peaks['dataset'] = 'Roadmap'\n",
    "df_immune_peaks.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_peaks = pd.concat([df_encode_peaks, df_immune_peaks]).reset_index(drop=True)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# Remove outliers, see 2023 - figures.ipynb\n",
    "OUTLIERS = [('H3K4me3', 'NK', ''),\n",
    "            ('H3K4me1', 'NK', ''),\n",
    "            ('H3K27me3', 'TCell', ''),\n",
    "            ('H3K27me3', 'NHLF', 'rep2'),\n",
    "            ('H3K27me3', 'GM12878', 'rep1')]\n",
    "df_peaks = df_peaks.loc[[(m, c, r) not in OUTLIERS\n",
    "                         for _, (m, c, r) in df_peaks[['modification', 'cell', 'replicate']].iterrows()]]"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(14, 4))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    sns.barplot(data=df_peaks[df_peaks['modification'] == m], x='tool', y='peaks', hue='dataset',\n",
    "                capsize=.2, err_kws={'linewidth': 2}, edgecolor=\"black\",\n",
    "                # palette=TOOLS_PALETTE,\n",
    "                order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                ax=ax)\n",
    "    sns.swarmplot(data=df_peaks[df_peaks['modification'] == m], x='tool', y='peaks',  hue='dataset',\n",
    "                  dodge=True, size=3, color=\"black\", alpha=0.5, palette='dark:black',\n",
    "                  order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                  legend=False,\n",
    "                  ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.title.set_text(m)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel(None)\n",
    "    else:\n",
    "        ax.set_ylabel('peaks')\n",
    "    if i < 4:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    # ax.set_ylim(0, 1.9e7)\n",
    "plt.tight_layout()\n",
    "# plt.savefig(GSE26320_PATH_HG38 + f'/pics/library_mapped_detailed_all.pdf', bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Compute BAM coverage"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def load_bams(path, modifications, cells, replicates, control='Input'):\n",
    "    df_bams = pd.DataFrame(columns=['bam_file', 'modification', 'cell', 'replicate', 'reads'], dtype=object)\n",
    "    for f in tqdm(os.listdir(path)):\n",
    "        if not f.endswith('.bam'):\n",
    "            continue\n",
    "        cell = next((c for c in cells if c in f), None)\n",
    "        rep = next((r for r in replicates if r in f), None)\n",
    "        mod = next((m for m in modifications if m in f), None)\n",
    "        if mod is None and control is not None and control in f:\n",
    "            mod = 'Control'\n",
    "        if cell and mod and rep is not None:\n",
    "            # # Limit cells\n",
    "            # if cell not in ['GM12878', 'K562', 'Huvec', 'TCell', 'BCell', 'Monocyte']:\n",
    "            #     continue\n",
    "            file = os.path.join(path, f)\n",
    "            reads = ! samtools view -F 0x04 -c {file}\n",
    "            reads = int(reads[0])\n",
    "            df_bams.loc[len(df_bams)] = (file, mod, cell, rep, reads)\n",
    "    return df_bams\n"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_encode_bams = load_bams(GSE26320_PATH + '/bams', MODIFICATIONS, GSE26320_CELLS, GSE26320_REPS, 'Input')\n",
    "df_encode_bams['dataset'] = 'ENCODE'\n",
    "df_encode_bams"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_immune_bams = load_bams(IMMUNE_PATH + '/bams', MODIFICATIONS, IMMUNE_CELLS, IMMUNE_REPS, 'Control')\n",
    "df_immune_bams['dataset'] = 'Roadmap'\n",
    "df_immune_bams"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_bams = pd.concat([df_encode_bams, df_immune_bams]).reset_index(drop=True)\n",
    "df_bams = df_bams.loc[[(m, c, r) not in OUTLIERS for _, (m, c, r) in df_bams[['modification', 'cell', 'replicate']].iterrows()]].copy()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(6, 4))\n",
    "ax = plt.axes()\n",
    "g_results = sns.barplot(data=df_bams,\n",
    "                        hue='modification', x='dataset', y='reads', ax=ax,\n",
    "                        capsize=.05, err_kws={'linewidth': 2}, edgecolor=\"black\",\n",
    "                        hue_order=MODIFICATIONS+['Control']\n",
    "                        )\n",
    "ax.set_title('Library size')\n",
    "# Put a legend to the right of the current axis\n",
    "g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Find control tracks"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "! mkdir -p {DOWNSAMPLE_DIR}/bams\n",
    "controls = []\n",
    "control_reads = []\n",
    "\n",
    "for _, (ds, m, c, r) in tqdm(\n",
    "        df_bams[['dataset', 'modification', 'cell', 'replicate']].iterrows()):\n",
    "    if m == 'Control':\n",
    "        controls.append(None)\n",
    "        control_reads.append(None)\n",
    "        continue\n",
    "    # # Limit cells\n",
    "    # if c not in ['GM12878', 'K562', 'Huvec', 'TCell', 'BCell', 'Monocyte']:\n",
    "    #     controls.append(None)\n",
    "    #     control_reads.append(None)\n",
    "    #     continue\n",
    "    print(ds, m, c, r)\n",
    "    t = df_bams[(df_bams['dataset'] == ds) & (df_bams['modification'] == 'Control') & (df_bams['cell'] == c)]\n",
    "    tr = t[t['replicate'] == r]\n",
    "    bam, reads = None, None\n",
    "    if len(tr) == 1:\n",
    "        bam = tr['bam_file'].values[0]\n",
    "        reads = tr['reads'].values[0]\n",
    "    elif len(t) == 1:\n",
    "        bam = t['bam_file'].values[0]\n",
    "        reads = t['reads'].values[0]\n",
    "    controls.append(bam)\n",
    "    control_reads.append(reads)\n",
    "    if bam is None:\n",
    "        print('No control', ds, m, c, r)\n",
    "        continue\n",
    "\n",
    "df_bams['control_file'] = controls\n",
    "df_bams['control_reads'] = control_reads\n",
    "df_bams.sample(3)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Mix control\n",
    "Use chromosome 1 for speed purposes."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "MIX_DIR = os.path.expanduser('~/data/2023_mix_chr1')\n",
    "MIX_SIGNAL_FRACTIONS = [0.7, 0.5, 0.2, 0.1]"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# ! mkdir -p {MIX_DIR}/bams\n",
    "# \n",
    "# tf = tempfile.mktemp()\n",
    "# \n",
    "# mixtures = []\n",
    "# for _, (ds, m, c, r, bam_file, control_file, reads, control_reads) in tqdm(\n",
    "#         df_bams[df_bams['modification'] != 'Control'][\n",
    "#             ['dataset', 'modification', 'cell', 'replicate', 'bam_file', 'control_file', 'reads', 'control_reads']\n",
    "#         ].iterrows()):\n",
    "#     # Limit cells\n",
    "#     if c not in ['GM12878', 'K562', 'Huvec', 'TCell', 'BCell', 'Monocyte']:\n",
    "#         continue    \n",
    "#     print(ds, m, c, r)\n",
    "#     control_name = os.path.basename(control_file)\n",
    "#     ! ln -sf {control_file} {MIX_DIR}/bams/{control_name}\n",
    "#     \n",
    "#     for fraction in MIX_SIGNAL_FRACTIONS:\n",
    "#         print(fraction)\n",
    "#         result_file = f'{MIX_DIR}/bams/{m}_{c}_{r}_{fraction}.bam'\n",
    "#         print(result_file)\n",
    "#         mixtures.append((ds, m, c, r, fraction, result_file))\n",
    "#         if os.path.exists(result_file):\n",
    "#             continue\n",
    "#         control_fraction = min(1.0, (1.0 - fraction) * reads / control_reads)\n",
    "#         ! samtools view -@ 4 -H {bam_file} > {tf}.sam\n",
    "#         ! samtools view -@ 4 -s {fraction} {bam_file} >> {tf}.sam\n",
    "#         ! samtools view -@ 4 -s {control_fraction} {control_file} >> {tf}.sam\n",
    "#         ! samtools view -S -b {tf}.sam > {tf}.bam\n",
    "#         ! samtools sort {tf}.bam -o {result_file}\n",
    "#         ! rm {tf}.sam {tf}.bam\n",
    "#     \n",
    "# df_mixtures = pd.DataFrame(mixtures, columns=['dataset', 'modification', 'cell', 'replicate', 'fraction', 'bam_file'])\n",
    "# del mixtures\n",
    "# df_mixtures.sample(3)        "
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "! mkdir -p {MIX_DIR}/bams\n",
    "\n",
    "tf = tempfile.mktemp()\n",
    "tab = '\\t'\n",
    "\n",
    "mixtures = []\n",
    "for _, (ds, m, c, r, bam_file, control_file, reads, control_reads) in tqdm(list(\n",
    "        df_bams[df_bams['modification'] != 'Control'][\n",
    "            ['dataset', 'modification', 'cell', 'replicate', 'bam_file', 'control_file', 'reads', 'control_reads']\n",
    "        ].iterrows())):\n",
    "    # # Limit cells\n",
    "    # if c not in ['GM12878', 'K562', 'Huvec', 'TCell', 'BCell', 'Monocyte']:\n",
    "    #     continue\n",
    "    print(ds, m, c, r)\n",
    "    control_name = os.path.basename(control_file)\n",
    "    if not os.path.exists(f'{MIX_DIR}/bams/{control_name}'):\n",
    "        print('Control')\n",
    "        ! samtools view -@ 4 -H {control_file} > {tf}.sam\n",
    "        ! samtools view -@ 4 {control_file} | grep \"chr1{tab}\" >> {tf}.sam\n",
    "        ! samtools view -S -b {tf}.sam > {tf}.bam\n",
    "        ! samtools sort {tf}.bam -o {MIX_DIR}/bams/{control_name}\n",
    "        ! rm {tf}.sam {tf}.bam\n",
    "        # ! ln -sf {control_file} {MIX_DIR}/bams/{control_name}\n",
    "\n",
    "    for fraction in MIX_SIGNAL_FRACTIONS:\n",
    "        print(fraction)\n",
    "        result_file = f'{MIX_DIR}/bams/{m}_{c}_{r}_{fraction}.bam'\n",
    "        print(result_file)\n",
    "        mixtures.append((ds, m, c, r, fraction, result_file))\n",
    "        if os.path.exists(result_file):\n",
    "            continue\n",
    "        control_fraction = min(1.0, (1.0 - fraction) * reads / control_reads)\n",
    "        ! samtools view -@ 4 -H {bam_file} > {tf}.sam\n",
    "        ! samtools view -@ 4 -s {fraction} {bam_file} | grep \"chr1{tab}\" >> {tf}.sam\n",
    "        ! samtools view -@ 4 -s {control_fraction} {control_file} | grep \"chr1{tab}\" >> {tf}.sam\n",
    "        ! samtools view -S -b {tf}.sam > {tf}.bam\n",
    "        ! samtools sort {tf}.bam -o {result_file}\n",
    "        ! rm {tf}.sam {tf}.bam\n",
    "\n",
    "df_mixtures = pd.DataFrame(mixtures, columns=['dataset', 'modification', 'cell', 'replicate', 'fraction', 'bam_file'])\n",
    "del mixtures\n",
    "df_mixtures.sample(3)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Peak calling and visualization\n",
    "Peak calling SPAN: important to keep fixed SPAN settings.\n",
    "\n",
    "```\n",
    "cd ~/data/2023_mi\n",
    "conda activate snakemake\n",
    "snakemake -p -s ~/work/chipseq-smk-pipeline/Snakefile all --cores all  --use-conda  --directory $(pwd) --config genome=hg38 fastq_dir=$(pwd) start_with_bams=True bams_dir=bams macs2=True span=False sicer=False --rerun-trigger mtime;\n",
    "snakemake -p -s ~/work/chipseq-smk-pipeline/Snakefile all --cores all  --use-conda  --directory $(pwd) --config genome=hg38 fastq_dir=$(pwd) start_with_bams=True bams_dir=bams macs2=True macs2_mode=broad macs2_params=\"--broad --broad-cutoff=0.1\" macs2_suffix=\"broad0.1\"  span=False sicer=False --rerun-trigger mtime; \n",
    "snakemake -p -s ~/work/chipseq-smk-pipeline/Snakefile all --cores all  --use-conda  --directory $(pwd) --config genome=hg38 fastq_dir=$(pwd) start_with_bams=True bams_dir=bams bamCoverage_params=\"--ignoreDuplicates\" span=True sicer=True span_bin=100 span_params=\"--keep-cache --noclip --debug  --blacklist /home/jetbrains/data/2023_wiggle/hg38-blacklist.v2.bed\" --rerun-trigger mtime\n",
    "```"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Analyze results"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "! mkdir -p {MIX_DIR}/original_chr1\n",
    "\n",
    "mixed_peaks_df = load_peaks_path(MIX_DIR, MODIFICATIONS, GSE26320_CELLS + IMMUNE_CELLS, GSE26320_REPS + IMMUNE_REPS)\n",
    "mixed_peaks_df['dataset'] = 'ENCODE'\n",
    "mixed_peaks_df.loc[mixed_peaks_df['cell'].isin(IMMUNE_CELLS), 'dataset'] = 'Roadmap'\n",
    "mixed_peaks_df['fraction'] = 0\n",
    "for fraction in MIX_SIGNAL_FRACTIONS:\n",
    "    mixed_peaks_df.loc[(mixed_peaks_df['file'].str.contains(f'_{fraction}_')) |\n",
    "                       (mixed_peaks_df['file'].str.contains(f'_{fraction}-')), 'fraction'] = fraction\n",
    "t = df_peaks[(df_peaks['cell'].isin(mixed_peaks_df['cell'].unique())) &\n",
    "             (df_peaks['replicate'].isin(mixed_peaks_df['replicate'].unique()))][[\n",
    "    'file', 'modification', 'cell', 'replicate', 'peaks', 'tool', 'dataset']].copy()\n",
    "t['fraction'] = 1.0\n",
    "# Processing original files\n",
    "t = df_peaks[(df_peaks['cell'].isin(mixed_peaks_df['cell'].unique())) &\n",
    "             (df_peaks['replicate'].isin(mixed_peaks_df['replicate'].unique()))][[\n",
    "    'file', 'modification', 'cell', 'replicate', 'peaks', 'tool', 'dataset']].copy()\n",
    "t['fraction'] = 1.0\n",
    "chr1_files = []\n",
    "chr1_peaks = []\n",
    "for file in tqdm(t['file']):\n",
    "    chr1_file = f'{MIX_DIR}/original_chr1/{os.path.basename(file)}'\n",
    "    if not os.path.exists(chr1_file):\n",
    "        tab = '\\t'\n",
    "        ! cat {file} | grep \"chr1{tab}\" > {chr1_file}\n",
    "    chr1_files.append(chr1_file)\n",
    "    chr1_peaks.append(lines(chr1_file))\n",
    "t['file'] = chr1_files\n",
    "t['peaks'] = chr1_peaks\n",
    "mixed_peaks_df = pd.concat([mixed_peaks_df, t]).reset_index(drop=True)\n",
    "mixed_peaks_df.sample(3)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Peaks number"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def plot_modified(df, y, line):\n",
    "    for ds in ['ENCODE', 'Roadmap']:\n",
    "        print(ds)\n",
    "        plt.figure(figsize=(15, 3))\n",
    "        axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "\n",
    "        for i, m in enumerate(MODIFICATIONS):\n",
    "            ax = axs[i]\n",
    "            t = df[(df['dataset'] == ds) & (df['modification'] == m)].copy()\n",
    "            fractions = sorted(df['fraction'].unique(), reverse=True)\n",
    "            t['j'] = 0\n",
    "            for j, fraction in enumerate(fractions):\n",
    "                t.loc[t['fraction'] == fraction, 'j'] = j\n",
    "            if line:\n",
    "                sns.lineplot(data=t,\n",
    "                             x='j', y=y, hue='tool',\n",
    "                             hue_order=['MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                             palette=TOOLS_PALETTE,\n",
    "                             ax=ax,\n",
    "                             errorbar='se', marker='o', markersize=10\n",
    "                             )\n",
    "                ax.set_xticks(range(len(fractions)), minor=False)\n",
    "                ax.set_xticklabels([str(f) for f in fractions])\n",
    "                ax.set_xlabel('fraction')\n",
    "                # We could use just invert axis, but this won't make linear dots\n",
    "                # ax.invert_xaxis()\n",
    "            else:\n",
    "                sns.barplot(data=t,\n",
    "                                    x='tool', y=y, hue='fraction', ax=ax,\n",
    "                                    capsize=.1, err_kws={'linewidth': 2}, edgecolor=\"black\",\n",
    "                                    order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                                    hue_order = sorted(df['fraction'].unique(), reverse=True),\n",
    "                                    )\n",
    "            ax.xaxis.set_tick_params(rotation=90)\n",
    "            ax.title.set_text(m)\n",
    "            if i > 0:\n",
    "                ax.set_ylabel(None)\n",
    "            if i < len(axs) - 1:\n",
    "                ax.legend().set_visible(False)\n",
    "            else:\n",
    "                # Put a legend to the right of the current axis\n",
    "                ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))    \n",
    "\n",
    "        plt.tight_layout()\n",
    "        plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Mixed')\n",
    "plot_modified(mixed_peaks_df, 'peaks', False)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "# mixed_peaks_df[mixed_peaks_df['tool'] == 'SPAN']",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Jaccard vs original peaks"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "import pyranges as pr\n",
    "\n",
    "def overlap_vs_original(df):\n",
    "    overlaps = []\n",
    "    for _, (ds, m, c, r, fraction, tool, file) in tqdm(list(df[df['fraction'] != 1.0][[\n",
    "        'dataset', 'modification', 'cell', 'replicate', 'fraction', 'tool', 'file'\n",
    "    ]].iterrows())):\n",
    "        t = df[(df['fraction'] == 1.0) &\n",
    "               (df['dataset'] == ds) &\n",
    "               (df['modification'] == m) &\n",
    "               (df['cell'] == c) &\n",
    "               (df['replicate'] == r) &\n",
    "               (df['tool'] == tool)]\n",
    "        if len(t) != 1:\n",
    "            print('No original peaks for', ds, m, c, r, tool)\n",
    "            continue\n",
    "        original_file = t['file'].values[0]\n",
    "        peaks = lines(file)\n",
    "        original_peaks = lines(original_file)\n",
    "        if peaks == 0 or original_peaks == 0:\n",
    "            peaks_len = int(bedl(file).sum())\n",
    "            original_peaks_len = int(bedl(original_file).sum())\n",
    "            overlaps.append((ds, m, c, r, fraction, tool,\n",
    "                             peaks, peaks_len, original_peaks, original_peaks_len,\n",
    "                             0, 0, 0, 0))\n",
    "            continue\n",
    "        file = pr.read_bed(file)\n",
    "        original_file = pr.read_bed(original_file)\n",
    "        original_peaks_len = original_file.lengths().sum()\n",
    "        peaks_len = file.lengths().sum()\n",
    "        # !bedtools intersect -a {original_file} -b {file} -wa -u > {tf}\n",
    "        original_peaks_overlap = len(original_file.overlap(file))\n",
    "        # !bedtools intersect -b {original_file} -a {file} -wa -u > {tf}\n",
    "        peaks_overlap = len(file.overlap(original_file))\n",
    "        # !bedtools intersect -a {original_file} -b {file} -wo > {tf}\n",
    "        overlap_len = file.intersect(original_file).lengths().sum()\n",
    "        overlaps.append((ds, m, c, r, fraction, tool, \n",
    "                         peaks, peaks_len, original_peaks, original_peaks_len,\n",
    "                         peaks_overlap, original_peaks_overlap, overlap_len, \n",
    "                         d(overlap_len, peaks_len + original_peaks_len - overlap_len)))\n",
    "        \n",
    "    return pd.DataFrame(overlaps, columns=[\n",
    "        'dataset', 'modification', 'cell', 'replicate', 'fraction', 'tool', \n",
    "        'peaks', 'peaks_len', 'original_peaks', 'original_peaks_len',\n",
    "        'peaks_overlap', 'original_peaks_overlap', 'overlap_len', 'jaccard'\n",
    "    ])"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "mixed_peaks_overlap_df = overlap_vs_original(mixed_peaks_df)\n",
    "mixed_peaks_overlap_df.sample(3)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "mixed_peaks_overlap_df[mixed_peaks_overlap_df['tool'] == 'SPAN']",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Add 1.0 and 0.0 fraction\n",
    "for _, (ds, m, c, r, t) in list(\n",
    "        mixed_peaks_overlap_df[\n",
    "            ['dataset', 'modification', 'cell', 'replicate', 'tool']].drop_duplicates().iterrows()):\n",
    "    mixed_peaks_overlap_df.loc[len(mixed_peaks_overlap_df)] = \\\n",
    "        [ds, m, c, r, 0.0, t] + [None] * (len(mixed_peaks_overlap_df.columns) - 7) + [0.0]\n",
    "    mixed_peaks_overlap_df.loc[len(mixed_peaks_overlap_df)] = \\\n",
    "        [ds, m, c, r, 1.0, t] + [None] * (len(mixed_peaks_overlap_df.columns) - 7) + [1.0]\n",
    "mixed_peaks_overlap_df.sample(3)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "from datetime import datetime\n",
    "now = datetime.now().isoformat(timespec='minutes')\n",
    "rfile = f'{MIX_DIR}/mixed_peaks_overlap_df_{now}.csv.gz'\n",
    "print(rfile)\n",
    "mixed_peaks_overlap_df.to_csv(rfile, compression='gzip', index=None)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Mixed control')\n",
    "plot_modified(mixed_peaks_overlap_df, 'jaccard', True)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('ENCODE')\n",
    "plt.figure(figsize=(4, 4))\n",
    "t = mixed_peaks_overlap_df[\n",
    "    (mixed_peaks_overlap_df['dataset']=='ENCODE') &\n",
    "    (mixed_peaks_overlap_df['modification']=='H3K4me1')\n",
    "    ].copy()\n",
    "t.sort_values(by=['fraction'], ascending=False, inplace=True)\n",
    "t['quality'] = [f'{int(f * 100)}%' for f in t['fraction']]\n",
    "sns.lineplot(data=t, x=\"quality\", y=\"jaccard\", hue='tool',\n",
    "             hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "             palette=TOOLS_PALETTE,\n",
    "             # errorbar=('sd', .5), \n",
    "             errorbar='se', \n",
    "             marker='o', markersize=8)\n",
    "plt.title('Jaccard vs original peaks')\n",
    "plt.xlabel('Quality')\n",
    "plt.ylabel('Jaccard')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Immune')\n",
    "plt.figure(figsize=(4, 4))\n",
    "ax = plt.axes()\n",
    "t = mixed_peaks_overlap_df[\n",
    "    (mixed_peaks_overlap_df['dataset']=='Roadmap') &\n",
    "    (mixed_peaks_overlap_df['modification']=='H3K4me1')\n",
    "    ].copy()\n",
    "t.sort_values(by=['fraction'], ascending=False, inplace=True)\n",
    "t['quality'] = [f'{int(f * 100)}%' for f in t['fraction']]\n",
    "sns.lineplot(data=t[t['tool'] != 'SPAN'], x=\"quality\", y=\"jaccard\", hue='tool',\n",
    "             hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "             palette=TOOLS_PALETTE,\n",
    "             errorbar='se',\n",
    "             marker='o', markersize=8,\n",
    "             ax=ax)\n",
    "sns.lineplot(data=t[t['tool'] == 'SPAN'], x=\"quality\", y=\"jaccard\", hue='tool',\n",
    "             hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "             palette=TOOLS_PALETTE,\n",
    "             errorbar='se',\n",
    "             marker='o', markersize=8,\n",
    "             legend=False,\n",
    "             ax=ax)\n",
    "\n",
    "plt.title('Jaccard vs original peaks')\n",
    "plt.xlabel('Quality')\n",
    "plt.ylabel('Jaccard')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "def compute_auc(df, x):\n",
    "    rows = []\n",
    "    for (ds, m, c, r, tool), t in df.groupby(['dataset', 'modification', 'cell', 'replicate', 'tool']):\n",
    "        t = t.copy()\n",
    "        t.sort_values(by=['fraction'], inplace=True)\n",
    "        auc = 0\n",
    "        fraction_prev = 0\n",
    "        value_prev = 0\n",
    "        for _, (fraction, value) in t[['fraction', x]].iterrows():\n",
    "            auc += (fraction - fraction_prev) * (value + value_prev) / 2\n",
    "            fraction_prev, value_prev = fraction, value\n",
    "        rows.append((ds, m, c, r, tool, auc))\n",
    "    return pd.DataFrame(data=rows, columns=['dataset', 'modification', 'cell', 'replicate', 'tool', 'auc'])"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "mixed_peaks_overlap_auc_df = compute_auc(mixed_peaks_overlap_df, 'jaccard')\n",
    "mixed_peaks_overlap_auc_df.sample(3)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "for ds in ['ENCODE', 'Roadmap']:\n",
    "    print(ds)\n",
    "    plt.figure(figsize=(6, 3))\n",
    "    ax = plt.axes()\n",
    "    ax.title.set_text('AUC jaccard vs original peaks')\n",
    "    sns.barplot(data=mixed_peaks_overlap_auc_df[mixed_peaks_overlap_auc_df['dataset']==ds], \n",
    "                x='modification', y='auc',\n",
    "                hue='tool',\n",
    "                capsize=.05, err_kws={'linewidth': 2},\n",
    "                order=MODIFICATIONS,\n",
    "                hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                ax=ax,\n",
    "                )\n",
    "    sns.swarmplot(data=mixed_peaks_overlap_auc_df[mixed_peaks_overlap_auc_df['dataset']==ds],\n",
    "                  x='modification', y='auc',\n",
    "                  hue='tool',\n",
    "                  dodge=True, size=2, color=\"black\", alpha=0.5, palette='dark:black',\n",
    "                  order=MODIFICATIONS,\n",
    "                  hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                  legend=False,\n",
    "                  ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=45)\n",
    "    # Put a legend to the right of the current axis\n",
    "    ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('ENCODE + Roadmap')\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "ax.title.set_text('AUC jaccard vs original peaks')\n",
    "sns.barplot(data=mixed_peaks_overlap_auc_df,\n",
    "            x='modification', y='auc',\n",
    "            hue='tool',\n",
    "            capsize=.2, err_kws={'linewidth': 2},\n",
    "            order=MODIFICATIONS,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            ax=ax,\n",
    "            )\n",
    "sns.swarmplot(data=mixed_peaks_overlap_auc_df,\n",
    "              x='modification', y='auc',\n",
    "              hue='tool',\n",
    "              dodge=True, size=1.5, color=\"black\", alpha=0.5, palette='dark:black',\n",
    "              order=MODIFICATIONS,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              legend=False,\n",
    "              ax=ax)\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "# Put a legend to the right of the current axis\n",
    "ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Jaccard between replicates"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def compute_reps_overlap(df_peaks, modifications, cells):\n",
    "    reps_overlap = pd.DataFrame(columns=['dataset', 'modification', 'cell', 'fraction', 'tool', 'rep1', 'rep2',\n",
    "                                         'peaks1', 'peaks1_len', 'peaks2', 'peaks2_len',\n",
    "                                         'peaks1_overlap', 'peaks2_overlap', 'peaks_overlap_len'], dtype=object)\n",
    "\n",
    "    tools = list(sorted(set(df_peaks['tool'])))\n",
    "    for ds, m, c, f in tqdm(list(product(['ENCODE', 'Roadmap'], modifications, cells, MIX_SIGNAL_FRACTIONS + [1.0]))):\n",
    "        print(ds, m, c, f)\n",
    "        tm = df_peaks[(df_peaks['dataset'] == ds) & (df_peaks['modification'] == m) & (df_peaks['cell'] == c) & (df_peaks['fraction'] == f)]\n",
    "        reps = list(sorted(set(tm['replicate'])))\n",
    "        for tool in tools:\n",
    "            # print(m, tool, len(reps))\n",
    "            files = [None] * len(reps)\n",
    "            peaks = [-1] * len(reps)\n",
    "            peaks_lens = [-1] * len(reps)\n",
    "            for i in range(len(reps)):\n",
    "                t1 = tm[(tm['tool'] == tool) & (tm['replicate'] == reps[i])]\n",
    "                if len(t1) == 0:\n",
    "                    continue\n",
    "                file = t1['file'].values[0]\n",
    "                peaks[i] = lines(file)\n",
    "                if peaks[i] == 0:\n",
    "                    peaks_lens[i] = 0\n",
    "                    continue\n",
    "                files[i] = pr.read_bed(file)\n",
    "                peaks_lens[i] = files[i].lengths().sum()\n",
    "            for i, j in product(range(len(reps)), range(len(reps))):\n",
    "                if i >= j:\n",
    "                    continue\n",
    "                rep1, rep2 = reps[i], reps[j]\n",
    "                file1, file2 = files[i], files[j]\n",
    "                peaks1, peaks2 = peaks[i], peaks[j]\n",
    "                peaks1_len, peaks2_len = peaks_lens[i], peaks_lens[j]\n",
    "                if peaks1 <= 0 or peaks2 <= 0:\n",
    "                    reps_overlap.loc[len(reps_overlap)] = \\\n",
    "                        (ds, m, c, f, tool, rep1, rep2, peaks1, peaks1_len, peaks2, peaks2_len,\n",
    "                         0, 0, 0)\n",
    "                    continue\n",
    "                overlap1 = len(file1.overlap(file2))\n",
    "                overlap2 = len(file2.overlap(file1))\n",
    "                overlap_len = file1.intersect(file2).lengths().sum()\n",
    "                reps_overlap.loc[len(reps_overlap)] = \\\n",
    "                    (ds, m, c, f, tool, rep1, rep2, peaks1, peaks1_len, peaks2, peaks2_len, \n",
    "                     overlap1, overlap2, overlap_len)\n",
    "    reps_overlap['jaccard'] = [\n",
    "        d(lo, l1 + l2 - lo)\n",
    "        for l1, l2, lo in zip(reps_overlap['peaks1_len'], reps_overlap['peaks2_len'], reps_overlap['peaks_overlap_len'])\n",
    "    ]\n",
    "    return reps_overlap\n"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_reps_overlap = compute_reps_overlap(mixed_peaks_df, MODIFICATIONS, GSE26320_CELLS + IMMUNE_CELLS)\n",
    "# Add 0.0 fraction\n",
    "for _, (ds, m, c, t) in list(\n",
    "        df_reps_overlap[\n",
    "            ['dataset', 'modification', 'cell', 'tool']].drop_duplicates().iterrows()):\n",
    "    df_reps_overlap.loc[len(df_reps_overlap)] = \\\n",
    "        [ds, m, c, 0.0, t] + [None] * (len(df_reps_overlap.columns) - 6) + [0.0]\n",
    "df_reps_overlap.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "mixed_peaks_df[(mixed_peaks_df['dataset'] == 'Roadmap') & \n",
    "               (mixed_peaks_df['modification'] == 'H3K27ac')][['cell', 'replicate']].drop_duplicates()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Mixed control replicates overlap')\n",
    "plot_modified(df_reps_overlap, 'jaccard', True)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('ENCODE')\n",
    "plt.figure(figsize=(4, 4))\n",
    "t = df_reps_overlap[\n",
    "    (df_reps_overlap['dataset']=='ENCODE') &\n",
    "    (df_reps_overlap['modification']=='H3K27ac') \n",
    "    ].copy()\n",
    "t.sort_values(by=['fraction'], ascending=False, inplace=True)\n",
    "t['quality'] = [f'{int(f * 100)}%' for f in t['fraction']]\n",
    "sns.lineplot(data=t, x=\"quality\", y=\"jaccard\", hue='tool',\n",
    "             hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "             palette=TOOLS_PALETTE,\n",
    "             # errorbar=('sd', .5),\n",
    "             errorbar='se', \n",
    "             marker='o', markersize=8)\n",
    "plt.title('Jaccard between replicates')\n",
    "plt.xlabel('Quality')\n",
    "plt.ylabel('Jaccard')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Immune')\n",
    "plt.figure(figsize=(4, 4))\n",
    "t = df_reps_overlap[\n",
    "    (df_reps_overlap['dataset']=='Roadmap') &\n",
    "    (df_reps_overlap['modification']=='H3K4me1')\n",
    "    ].copy()\n",
    "t.sort_values(by=['fraction'], ascending=False, inplace=True)\n",
    "t['quality'] = [f'{int(f * 100)}%' for f in t['fraction']]\n",
    "sns.lineplot(data=t, x=\"quality\", y=\"jaccard\", hue='tool',\n",
    "             hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "             palette=TOOLS_PALETTE,\n",
    "             # errorbar=('sd', .5),\n",
    "             errorbar='se', \n",
    "             marker='o', markersize=8)\n",
    "plt.title('Jaccard between replicates')\n",
    "plt.xlabel('Quality')\n",
    "plt.ylabel('Jaccard')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "def compute_auc_reps(df, x):\n",
    "    rows = []\n",
    "    for (ds, m, c, tool), t in df.groupby(['dataset', 'modification', 'cell', 'tool']):\n",
    "        t = t.copy()\n",
    "        t.sort_values(by=['fraction'], inplace=True)\n",
    "        auc = 0\n",
    "        fraction_prev = 0\n",
    "        value_prev = 0\n",
    "        for _, (fraction, value) in t[['fraction', x]].iterrows():\n",
    "            auc += (fraction - fraction_prev) * (value + value_prev) / 2\n",
    "            fraction_prev, value_prev = fraction, value\n",
    "        rows.append((ds, m, c, tool, auc))\n",
    "    return pd.DataFrame(data=rows, columns=['dataset', 'modification', 'cell', 'tool', 'auc'])"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_reps_overlap_auc_df = compute_auc_reps(df_reps_overlap, 'jaccard')\n",
    "df_reps_overlap_auc_df.sample(3)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "for ds in ['ENCODE', 'Roadmap']:\n",
    "    print(ds)\n",
    "    plt.figure(figsize=(6, 3))\n",
    "    ax = plt.axes()\n",
    "    ax.title.set_text('AUC jaccard between replicates')\n",
    "    sns.barplot(data=df_reps_overlap_auc_df[df_reps_overlap_auc_df['dataset']==ds],\n",
    "                x='modification', y='auc',\n",
    "                hue='tool',\n",
    "                capsize=.05, err_kws={'linewidth': 2},\n",
    "                order=MODIFICATIONS,\n",
    "                hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                ax=ax,\n",
    "                )\n",
    "    sns.swarmplot(data=df_reps_overlap_auc_df[df_reps_overlap_auc_df['dataset']==ds],\n",
    "                  x='modification', y='auc',\n",
    "                  hue='tool',\n",
    "                  dodge=True, size=2, color=\"black\", alpha=0.5, palette='dark:black',\n",
    "                  order=MODIFICATIONS,\n",
    "                  hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                  legend=False,\n",
    "                  ax=ax)\n",
    "\n",
    "    ax.xaxis.set_tick_params(rotation=45)\n",
    "    # Put a legend to the right of the current axis\n",
    "    ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('ENCODE + Roadmap')\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "ax.title.set_text('AUC jaccard between replicates')\n",
    "sns.barplot(data=df_reps_overlap_auc_df,\n",
    "            x='modification', y='auc',\n",
    "            hue='tool',\n",
    "            capsize=.2, err_kws={'linewidth': 2},\n",
    "            order=MODIFICATIONS,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            ax=ax,\n",
    "            )\n",
    "sns.swarmplot(data=df_reps_overlap_auc_df,\n",
    "              x='modification', y='auc',\n",
    "              hue='tool',\n",
    "              dodge=True, size=2, color=\"black\", alpha=0.5, palette='dark:black',\n",
    "              order=MODIFICATIONS,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              legend=False,\n",
    "              ax=ax)\n",
    "\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "# Put a legend to the right of the current axis\n",
    "ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}