2018 - ABF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2018 ABF SPAN\n",
    "\n"
   ]
  },
  {
   "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 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 tempfile"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "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"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def load_peaks(path, suffix, modifications, replicates):\n",
    "    df = pd.DataFrame(columns=['file', 'modification', 'replicate', 'peaks'],\n",
    "                      dtype=object)\n",
    "    for f in tqdm(os.listdir(path)):\n",
    "        if not f.endswith(suffix):\n",
    "            continue\n",
    "        mod = next((m for m in modifications if m.lower() in f.lower()), None)\n",
    "        rep = next((r for r in replicates if f'{r.lower()}_' in f.lower()), None)\n",
    "        if 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, rep, ps)\n",
    "    return df\n"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Load peaks"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "PATH = os.path.expanduser('~/data/2018_chipseq_y20o20')\n",
    "MODIFICATIONS = ['k4me3', 'k27ac', 'k4me1', 'k27me3', 'k36me3']\n",
    "REPLICATES = [f'OD{i}' for i in range(1, 30)] + [f'YD{i}' for i in range(1, 30)]"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_macs2 = load_peaks(os.path.join(PATH, 'macs2'), '.narrowPeak', MODIFICATIONS, 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, 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'), '-FDR0.01', MODIFICATIONS, 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, REPLICATES, )\n",
    "df_span['tool'] = 'SPAN'\n",
    "print('SPAN', len(df_span))\n",
    "\n",
    "df_macs2broad_abf = load_peaks(os.path.join(PATH, 'macs2_abf'), '.broadPeak', MODIFICATIONS, REPLICATES, )\n",
    "df_macs2broad_abf['tool'] = 'MACS2 ABF'\n",
    "print('MACS2 broad ABF', len(df_macs2broad_abf))\n",
    "\n",
    "df_sicer_abf = load_peaks(os.path.join(PATH, 'sicer_abf'), '', MODIFICATIONS, REPLICATES, )\n",
    "df_sicer_abf['tool'] = 'SICER ABF'\n",
    "print('SICER ABF', len(df_sicer_abf))\n",
    "\n",
    "df_span_tuned = load_peaks(os.path.join(PATH, 'span_tuned'), '.bed', MODIFICATIONS, REPLICATES, )\n",
    "df_span_tuned['tool'] = 'SPAN tuned'\n",
    "print('SPAN tuned', len(df_span_tuned))\n",
    "\n",
    "df_peaks = pd.concat([df_macs2broad_abf, df_sicer_abf, df_span_tuned,\n",
    "                      df_macs2, df_macs2broad,  df_sicer, df_span]).reset_index(drop=True)\n",
    "df_peaks.sample(3)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(12, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    dfm = df_peaks[df_peaks['modification'] == m]\n",
    "    sns.barplot(data=dfm, x='tool', y='peaks',\n",
    "                order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                capsize=.2, err_kws={'linewidth': 2}, edgecolor=\"black\",\n",
    "                ax=ax)\n",
    "    sns.stripplot(data=dfm, x='tool', y='peaks',\n",
    "                  order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                  size=1, color='black', alpha=0.5,\n",
    "                  ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.title.set_text(m)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# Drop outliers\n",
    "OUTLIERS = {'k27ac': ['YD1', 'YD6'],\n",
    "           'k27me3': ['OD9', 'YD1', 'YD2', 'YD3', 'YD4', 'YD9', 'YD10', 'YD11'],\n",
    "           'k36me3': ['OD3', 'OD6', 'OD12', 'OD18', 'OD20', 'YD1', 'YD3', 'YD4', 'YD5'],\n",
    "           'k4me3': ['OD6', 'OD7', 'OD14', 'YD2', 'YD3', 'YD10', 'YD14']}\n",
    "t_outliers = [any(f'{o}_h3{m}'.lower() in f.lower() or f'{o}_{m}'.lower() in f.lower() \n",
    "                  for m, ol in OUTLIERS.items() for o in ol)\n",
    "              for f in df_peaks['file']]\n",
    "print('Total outliers', sum(t_outliers))\n",
    "df_peaks_no = df_peaks[[not to for to in t_outliers]].copy()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "df_peaks_no['tool'].unique()",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# plt.figure(figsize=(1, 2))\n",
    "# sns.boxplot(data=df_peaks_no, x='tool', y='peaks',\n",
    "#             # order=['MACS2', 'SICER', 'SPAN tuned', 'SPAN 0.05', 'SPAN 1e-4'],\n",
    "#             )\n",
    "# # plt.title('K36me3')\n",
    "# plt.tight_layout()\n",
    "# plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(12, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    dfm = df_peaks_no[df_peaks_no['modification'] == m]\n",
    "    # sns.boxplot(data=dfm, x='tool', y='peaks',\n",
    "    #             order=['MACS2', 'SICER', 'SPAN tuned', 'SPAN'],\n",
    "    #             showfliers=True,\n",
    "    #             ax=ax)\n",
    "    sns.barplot(data=dfm, x='tool', y='peaks',\n",
    "                order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                capsize=.2, edgecolor=\"black\",\n",
    "                err_kws={'linewidth': 2},\n",
    "                ax=ax)\n",
    "    sns.stripplot(data=dfm, x='tool', y='peaks',\n",
    "                  order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                  size=1.5, alpha=0.5, color='black',\n",
    "                  ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.title.set_text(m)\n",
    "    ax.set_ylim(0, min(1.2e5, dfm['peaks'].max()))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Lengths"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('Load lengths')\n",
    "ts = []\n",
    "for mod, rep, file, tool in tqdm(zip(\n",
    "        df_peaks_no['modification'], df_peaks_no['replicate'], df_peaks_no['file'],\n",
    "        df_peaks_no['tool']\n",
    ")):\n",
    "    lengths = bedl(file)\n",
    "    t = pd.DataFrame(dict(length=lengths))\n",
    "    t = t.sample(min(len(t), 10_000)).copy()\n",
    "    t['modification'] = mod\n",
    "    t['replicate'] = rep\n",
    "    t['tool'] = tool\n",
    "    t['footprint'] = sum(lengths)\n",
    "    ts.append(t)\n",
    "df_lens = pd.concat(ts).reset_index(drop=True)\n",
    "del ts, t\n",
    "df_lens.sample(10)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sns.set_theme(style=\"whitegrid\")\n",
    "plt.figure(figsize=(12, 3.5))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    sns.boxplot(data=df_lens[df_lens['modification'] == m], x='tool', y='length',\n",
    "                showfliers=False,\n",
    "                order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'SICER', 'SPAN'],\n",
    "                ax=ax)\n",
    "    ax.title.set_text(m)\n",
    "    if m not in ['H3K36me3']:\n",
    "        ax.set_ylim(bottom=20, top=1e5)\n",
    "    else:\n",
    "        ax.set_ylim(bottom=20, top=1e5)\n",
    "    ax.set(yscale='log')\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Jaccard"
  },
  {
   "metadata": {
    "collapsed": false
   },
   "cell_type": "code",
   "source": [
    "from itertools import product\n",
    "import pyranges as pr\n",
    "\n",
    "def compute_reps_overlap(df_peaks, modifications, fixed=-1):\n",
    "    reps_overlap = pd.DataFrame(columns=['modification', '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 m in modifications:\n",
    "        tm = df_peaks[(df_peaks['modification'] == m)]\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",
    "            # print('Load')\n",
    "            for i in range(len(reps)):\n",
    "                t1 = tm[(tm['tool'] == tool) & (tm['replicate'] == reps[i])]\n",
    "                if len(t1) == 0:\n",
    "                    continue\n",
    "                file = t1['file'].values[0]\n",
    "                peaks[i] = lines(file)\n",
    "                if peaks[i] == 0:\n",
    "                    peaks_lens[i] = 0\n",
    "                    continue\n",
    "                files[i] = pr.read_bed(sorted_file(file))\n",
    "                peaks_lens[i] = files[i].lengths().sum()\n",
    "            # print('Intersect')\n",
    "            for i, j in tqdm(list(product(range(len(reps)), range(len(reps))))):\n",
    "                if fixed != -1 and i != fixed or i >= j:\n",
    "                    continue\n",
    "                rep1, rep2 = reps[i], reps[j]\n",
    "                file1, file2 = files[i], files[j]\n",
    "                peaks1, peaks2 = peaks[i], peaks[j]\n",
    "                peaks1_len, peaks2_len = peaks_lens[i], peaks_lens[j]\n",
    "                if peaks1 <= 0 or peaks2 <= 0:\n",
    "                    reps_overlap.loc[len(reps_overlap)] = \\\n",
    "                        (m, tool, rep1, rep2, peaks1, peaks1_len, peaks2, peaks2_len, \n",
    "                         0, 0, 0)\n",
    "                    continue\n",
    "                overlap1 = len(file1.overlap(file2))\n",
    "                overlap2 = len(file2.overlap(file1))\n",
    "                overlap_len = file1.intersect(file2).lengths().sum()\n",
    "                # assert file2.intersect(file1).lengths().sum() == overlap_len\n",
    "                reps_overlap.loc[len(reps_overlap)] = \\\n",
    "                    (m, 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"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "df_overlap_full = compute_reps_overlap(df_peaks_no, MODIFICATIONS)",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "sns.boxplot(data=df_overlap_full, x='modification', hue='tool', y='jaccard',\n",
    "            hue_order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "            showfliers=False,\n",
    "            ax=ax)\n",
    "ax.set_title('Jaccard between replicates')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "t = pd.concat([df_overlap_full[['modification', 'tool']],\n",
    "               df_overlap_full[['modification', 'tool']]]).reset_index(drop=True)\n",
    "overlaps = [d(o, p) for o, p in zip(df_overlap_full['peaks1_overlap'], df_overlap_full['peaks1'])] + \\\n",
    "           [d(o, p) for o, p in zip(df_overlap_full['peaks2_overlap'], df_overlap_full['peaks2'])]\n",
    "t['overlap'] = overlaps\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "sns.boxplot(data=t, x='modification', hue='tool', y='overlap',\n",
    "            hue_order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "            showfliers=False\n",
    "            )\n",
    "ax.set_title('Overlap between replicates')\n",
    "ax.set_ylabel('Overlap')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Comparison vs ENCODE CD14 data"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "```\n",
    "cd ~/data/2024_cd14_monocytes_chipseq\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 span=True span_params=\"--debug --keep-cache\" span_threads=2 --rerun-trigger mtime \n",
    "```\n",
    "and use liftover to hg19\n",
    "```\n",
    "for F in sicer/*hg38*FDR0.01; do echo $F; NF=${F/hg38/hg19}; echo $NF; liftOver -bedPlus=6 $F hg38ToHg19.over.chain $NF ${NF}_unmapped; done\n",
    "for F in macs2/*hg38*Peak; do echo $F; NF=${F/hg38/hg19}; echo $NF; liftOver -bedPlus=6 $F hg38ToHg19.over.chain $NF ${NF}_unmapped; done\n",
    "for F in span/*hg38*.peak; do echo $F; NF=${F/hg38/hg19}; echo $NF; liftOver -bedPlus=6 $F hg38ToHg19.over.chain $NF ${NF}_unmapped; done\n",
    "```"
   ]
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "CD14_ENCODE_PATH = os.path.expanduser('~/data/2024_cd14_monocytes_chipseq')\n",
    "df_cd14_macs2 = load_peaks(os.path.join(CD14_ENCODE_PATH, 'macs2'), '.narrowPeak', MODIFICATIONS, [''], )\n",
    "df_cd14_macs2['tool'] = 'MACS2'\n",
    "print('MACS2', len(df_cd14_macs2))\n",
    "\n",
    "df_cd14_macs2broad = load_peaks(os.path.join(CD14_ENCODE_PATH, 'macs2'), '.broadPeak', MODIFICATIONS, [''], )\n",
    "df_cd14_macs2broad['tool'] = 'MACS2 broad'\n",
    "print('MACS2 broad', len(df_cd14_macs2broad))\n",
    "\n",
    "df_cd14_sicer = load_peaks(os.path.join(CD14_ENCODE_PATH, 'sicer'), '-FDR0.01', MODIFICATIONS, [''], )\n",
    "df_cd14_sicer['tool'] = 'SICER'\n",
    "print('SICER', len(df_cd14_sicer))\n",
    "\n",
    "df_cd14_span = load_peaks(os.path.join(CD14_ENCODE_PATH, 'span'), '.peak', MODIFICATIONS, [''], )\n",
    "df_cd14_span['tool'] = 'SPAN'\n",
    "print('SPAN', len(df_cd14_span))\n",
    "\n",
    "df_cd14 = pd.concat([df_cd14_macs2, df_cd14_macs2broad, df_cd14_sicer, df_cd14_span]).reset_index(drop=True)\n",
    "# Ignore multiple versions for several modifications\n",
    "df_cd14 = df_cd14[(df_cd14['file'].str.contains('GSM1102807')) & (df_cd14['file'].str.contains('hg19'))]\n",
    "df_cd14.sample(3)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(10, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    dfm = df_cd14[df_cd14['modification'] == m]\n",
    "    sns.barplot(data=dfm, x='tool', y='peaks',\n",
    "                order=['MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                capsize=.2, edgecolor=\"black\",\n",
    "                err_kws={'linewidth': 2},\n",
    "                ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.title.set_text(m)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_peaks_no['dataset'] = 'Y20O20'\n",
    "df_cd14['dataset'] = 'ENCODE'\n",
    "t = pd.concat([df_peaks_no, df_cd14]).reset_index(drop=True)\n",
    "\n",
    "plt.figure(figsize=(14, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    dfm = t[t['modification'] == m]\n",
    "    sns.barplot(data=dfm, x='tool', y='peaks', hue='dataset',\n",
    "                order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                capsize=.2, edgecolor=\"black\",\n",
    "                err_kws={'linewidth': 2},\n",
    "                ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.title.set_text(m)\n",
    "    # ax.set_ylim(0, min(1.5e5, dfm['peaks'].max()))\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()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "t = df_cd14[df_cd14['tool'].isin(['MACS2', 'SPAN'])].copy()\n",
    "t['replicate'] = t['tool']\n",
    "t['tool'] = 'MACS2&SPAN'\n",
    "df_cd14_overlap = compute_reps_overlap(t, MODIFICATIONS)\n",
    "df_cd14_overlap"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('ENCODE')\n",
    "plt.figure(figsize=(5, 3))\n",
    "ax = plt.axes()\n",
    "sns.barplot(data=df_cd14_overlap, x='modification', hue='tool', y='jaccard',\n",
    "            ax=ax)\n",
    "ax.set_title('Jaccard between CD14 ENCODE MACS2 & SPAN')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "rows = []\n",
    "for _, row in df_cd14_overlap.iterrows():\n",
    "    rows.append((row['modification'], 'macs2_by_span', d(row['peaks1_overlap'], row['peaks1'])))\n",
    "    rows.append((row['modification'], 'span_by_macs2', d(row['peaks2_overlap'], row['peaks2'])))\n",
    "t = pd.DataFrame(rows, columns=['modification', 'name', 'overlap'])\n",
    "del rows\n",
    "plt.figure(figsize=(5, 3))\n",
    "ax = plt.axes()\n",
    "sns.barplot(data=t, x='modification', hue='name', y='overlap',\n",
    "            ax=ax)\n",
    "ax.set_title('Overlap between CD14 ENCODE MACS2 & SPAN')\n",
    "ax.set_ylabel('Overlap')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for (m, tool), dfmt in df_peaks_no.groupby(['modification', 'tool']):\n",
    "    if tool in ['MACS2', 'MACS2 ABF']:\n",
    "        encode_tool = 'MACS2'\n",
    "    elif tool == 'MACS2 broad':\n",
    "        encode_tool = 'MACS2 broad'\n",
    "    elif tool in ['SICER', 'SICER ABF']:\n",
    "        encode_tool = 'SICER'\n",
    "    elif tool in ['SPAN', 'SPAN tuned']:\n",
    "        encode_tool = 'SPAN'\n",
    "    else:\n",
    "        continue\n",
    "    t = df_cd14[(df_cd14['modification'] == m) & (df_cd14['tool'] == encode_tool)].copy().reset_index(drop=True)\n",
    "    t['replicate'] = 'ENCODE'\n",
    "    t['tool'] = tool\n",
    "    ts.append(t)\n",
    "    ts.append(dfmt.copy().reset_index(drop=True))\n",
    "df_joint = pd.concat(ts).reset_index(drop=True)\n",
    "df_joint.head(5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_joint_overlap = compute_reps_overlap(df_joint, MODIFICATIONS, fixed=0)\n",
    "df_joint_overlap.sample(3)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(8, 3))\n",
    "ax = plt.axes()\n",
    "sns.barplot(data=df_joint_overlap, x='modification', y='jaccard', hue='tool',\n",
    "            hue_order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "            ax=ax)\n",
    "ax.set_title('Jaccard between Y20O20 and CD14 ENCODE ')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "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'\n",
    "\n",
    "\n",
    "TOOLS = ['MACS2', 'MACS2 broad', 'SPAN', 'SICER']\n",
    "t = df_joint_overlap[df_joint_overlap['tool'].isin(TOOLS)].copy()\n",
    "\n",
    "palette = plt.get_cmap('tab10')\n",
    "TOOLS_PALETTE = {t: palette(i) for i, t in enumerate(TOOLS)}\n",
    "\n",
    "update_abf_modifications(t)\n",
    "\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "sns.barplot(data=t, x='modification', y='jaccard', hue='tool',\n",
    "            hue_order=TOOLS,\n",
    "            palette=TOOLS_PALETTE,\n",
    "            capsize=.2, err_kws={'linewidth': 2},\n",
    "            ax=ax)\n",
    "sns.stripplot(data=t, x='modification', y='jaccard', hue='tool',\n",
    "                  dodge=True, size=2, palette='dark:black', alpha=0.4,\n",
    "              legend=False,\n",
    "                  hue_order=TOOLS,\n",
    "                  ax=ax)\n",
    "ax.set_title('Jaccard between ULI and ENCODE')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_joint_overlap['precision'] = [d(po, p) for po, p in zip(df_joint_overlap['peaks1_overlap'], df_joint_overlap['peaks1'])]\n",
    "df_joint_overlap['recall'] = [d(po, p) for po, p in zip(df_joint_overlap['peaks2_overlap'], df_joint_overlap['peaks2'])]\n",
    "df_joint_overlap['f1'] = [2 / (1 / p + 1 / r) if min(p, r) > 0 else 0 for p, r in zip(df_joint_overlap['precision'], df_joint_overlap['recall'])]"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(8, 3))\n",
    "ax = plt.axes()\n",
    "sns.barplot(data=df_joint_overlap, x='modification', y='f1', hue='tool',\n",
    "            hue_order=['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "            ax=ax)\n",
    "ax.set_title('F1 between Y20O20 and CD14 ENCODE ')\n",
    "ax.set_ylabel('F1')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Join\n",
    "rows = []\n",
    "for _, row in df_joint_overlap.iterrows():\n",
    "    rows.append((row['modification'], row['tool'], 'ENCODE_by_Y20O20', d(row['peaks1_overlap'], row['peaks1'])))\n",
    "    rows.append((row['modification'], row['tool'], 'Y20O20_by_ENCODE', d(row['peaks2_overlap'], row['peaks2'])))\n",
    "overlap_df = pd.DataFrame(rows, columns=['modification', 'tool', 'name', 'overlap'])\n",
    "del rows\n",
    "overlap_df.sample(3)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(20, 3))\n",
    "axs = [plt.subplot(1, 7, i + 1) for i in range(7)]\n",
    "for i, tool in enumerate(['MACS2 ABF', 'SICER ABF', 'SPAN tuned', 'MACS2', 'MACS2 broad', 'SICER', 'SPAN']):\n",
    "    print(tool)\n",
    "    ax = axs[i]\n",
    "    sns.barplot(data=overlap_df[overlap_df['tool'] == tool], x='modification', hue='name', y='overlap',\n",
    "                ax=ax)\n",
    "    ax.set_title(f'Overlap {tool}')\n",
    "    ax.set_ylabel('Overlap')\n",
    "    ax.set_xlabel('Modification')\n",
    "    ax.set_ylim(0, 1)\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",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Consensus peaks vs ChromHMM markup"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "CONSENSUS_PEAKS = {\n",
    "    m: f'{PATH}/{m}_consensus_union50%.bed' for m in MODIFICATIONS\n",
    "}\n",
    "for m in MODIFICATIONS:\n",
    "    ! wc -l {CONSENSUS_PEAKS[m]}\n",
    "    pass"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('Split chromHMM into separate files')\n",
    "# CHROMHMM_FILE = f'{PATH}/cd14_chromhmm_hg19_ENCFF808WWL.bed'\n",
    "CHROMHMM_FILE = f'{PATH}/cd14_chromhmm_hg19_ENCFF072SCA.bed'\n",
    "CHROMM_STATES = {}\n",
    "t = pd.read_csv(CHROMHMM_FILE, sep='\\t', header=None, usecols=[3])[3]\n",
    "for state in t.unique():\n",
    "    state = state.replace('/', '_')\n",
    "    f = CHROMHMM_FILE.replace('.bed', f'_{state}.bed')\n",
    "    CHROMM_STATES[state] = f\n",
    "    ! cat {CHROMHMM_FILE} | grep {state} > {f}"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "tf = tempfile.mktemp()\n",
    "\n",
    "rows = []\n",
    "for m in ['k4me1', 'k4me3', 'k27ac', 'k27me3', 'k36me3']:\n",
    "    print(m)\n",
    "    peaks_file = sorted_file(CONSENSUS_PEAKS[m]) \n",
    "    peaks = lines(peaks_file)\n",
    "    row = [m, peaks]\n",
    "    for state, state_file in tqdm(CHROMM_STATES.items()):\n",
    "        ! bedtools intersect -b {peaks_file} -a {state_file} -wa -u > {tf}\n",
    "        state_peaks = lines(tf)\n",
    "        row.append(state_peaks)\n",
    "    rows.append(row)\n",
    "\n",
    "chromhmm_ovlp_df = pd.DataFrame(rows, columns=['modification', 'peaks'] + list(CHROMM_STATES.keys()))\n",
    "chromhmm_ovlp_df"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import re\n",
    "chromhmm_rel_df = chromhmm_ovlp_df.copy()\n",
    "for state in CHROMM_STATES.keys():\n",
    "    chromhmm_rel_df[state] = chromhmm_rel_df[state] / lines(CHROMM_STATES[state]) if lines(CHROMM_STATES[state]) > 0 else 0\n",
    "chromhmm_rel_df.set_index('modification', inplace=True)\n",
    "chromhmm_rel_df.rename({c: re.sub('[0-9]+_', '', c) for c in chromhmm_rel_df.columns}, axis=1, inplace=True)\n",
    "chromhmm_rel_df = chromhmm_rel_df[['TssA', 'TssFlnk', 'TssFlnkU', 'TssFlnkD', 'Tx', 'TxWk',\n",
    "                                   'EnhG1', 'EnhG2', 'EnhA1', 'EnhA2', 'EnhWk', 'ZNF_Rpts',\n",
    "                                   'Het', 'TssBiv', 'EnhBiv', 'ReprPC', 'ReprPCWk', 'Quies']]\n",
    "chromhmm_rel_df"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(6, 2.5))\n",
    "ax = plt.axes()\n",
    "sns.heatmap(chromhmm_rel_df, cmap='coolwarm', ax=ax)\n",
    "plt.title('Overlap of consensus peaks with ChromHMM')\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "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
}