SPAN Comparison Opendnase.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison of span vs other peak calling tools\n",
    "\n",
    "Method based on F score described in original paper:\n",
    "Koohy H, Down TA, Spivakov M, Hubbard T (2014) \n",
    "A Comparison of Peak Callers Used for DNase-Seq Data. \n",
    "PLoS ONE 9(5): e96303. doi:10.1371/journal.pone.0096303\n",
    "\n",
    "http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0096303\n",
    "\n",
    "## Requirements\n",
    "\n",
    "* span.jar\n",
    "* java8\n",
    "* bedtools\n",
    "* wget\n",
    "* awk\n",
    "* ~20G free space to download all the required bed and bam files"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm.auto import tqdm\n",
    "import os\n",
    "import sys\n",
    "from os import listdir\n",
    "def run(cmd):\n",
    "    print(cmd)\n",
    "    if os.system(cmd) != 0:\n",
    "        sys.exit(1)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download all the \"Reference peaks set\" mentioned in http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0096303\n",
    "\n"
   ]
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "PATH = os.path.expanduser('~/data/2023_opendnase')\n",
    "RESOURCE = 'http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhTfbs'\n",
    "CELLS = ['GM12878', 'HelaS3', 'K562']"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "! mkdir ~/data/2023_opendnase\n",
    "! cd ~/data/2023_opendnase && wget \"http://journals.plos.org/plosone/article/asset?unique&id=info:doi/10.1371/journal.pone.0096303.s007\" -O supplementary.zip\n",
    "! cd ~/data/2023_opendnase && unzip -x   supplementary.zip"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have Supplemetary_Data/GM12878.txt, etc, we can download all the peaks automatically.\n",
    "Files are available here: http://journals.plos.org/plosone/article/asset?unique&id=info:doi/10.1371/journal.pone.0096303.s007"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "import pandas as pd\n",
    "\n",
    "for cell in tqdm(CELLS):\n",
    "    print(cell)\n",
    "    os.mkdir(PATH + '/' + cell)\n",
    "    df = pd.read_csv(PATH + f'/Supplementary_Data/{cell}.txt', sep=' ', header=None)\n",
    "    df.columns = ['file', 'checksum']\n",
    "    for f in df['file']:\n",
    "        ! cd {PATH} && wget {RESOURCE}/{f} -P {cell}"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we should all the reference peaks using Bedops as described in the article, however Bedops FAILS with segfault.\n",
    "\n",
    "Bedtools for the rescue!"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Lets merge them and filter chr22 to get reference peaks\n",
    "CELLS = ['GM12878', 'HelaS3', 'K562']\n",
    "! mkdir {PATH}/peaks\n",
    "for cell in tqdm(CELLS):\n",
    "    ! cd {PATH} && cat {cell}/*.gz > {cell}_all.bed\n",
    "    ! bedtools sort -i {PATH}/{cell}_all.bed > {PATH}/{cell}_sorted.bed\n",
    "    ! bedtools merge -i {PATH}/{cell}_sorted.bed > {PATH}/{cell}_peaks.bed\n",
    "    run(f'awk \\'$1==\"chr22\" {{print}}\\' {PATH}/{cell}_peaks.bed > {PATH}/peaks/{cell}_peaks_chr22.bed')\n",
    "! rm *_all.bed *_sorted.bed *_peaks.bed"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download and convert reads."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "! wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwDnase/wgEncodeUwDnaseK562AlnRep1.bam -O {PATH}/K562.bam\n",
    "! wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwDnase/wgEncodeUwDnaseGm12878AlnRep1.bam -O {PATH}/GM12878.bam\n",
    "! wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwDnase/wgEncodeUwDnaseHelas3AlnRep1.bam -O {PATH}/HelaS3.bam\n",
    "\n",
    "print('Lets convert downloaded bam files to bed')\n",
    "for bam in tqdm([f for f in listdir(PATH) if f.endswith(\".bam\")]):\n",
    "    ! bamToBed -bed12 -i {PATH}/{bam} > {PATH}/{bam}.bed"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can launch span to get enrichment for obtained bed tracks."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "! cd {PATH} && wget https://download.jetbrains.com/biolabs/span/span-2.0.6610.jar -O span.jar\n",
    "! cd {PATH} && wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.chrom.sizes"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "! mkdir {PATH}/span\n",
    "FDRS = [1e-40, 1e-20, 1e-15, 1e-10, 1e-8, 1e-6, 1e-4, 1e-3, 0.05, 0.1, 0.2, 0.5]\n",
    "for reads in tqdm([f for f in listdir(PATH) if f.endswith(\".bam.bed\")]):\n",
    "    for fdr in FDRS:\n",
    "        print(f'Processing {reads} FDR={fdr}')\n",
    "        ! cd {PATH} && java -jar span.jar analyze -cs hg19.chrom.sizes -t {reads} --peaks span/{reads}_{fdr}_span.bed --bin 100 --fdr {fdr} --keep-cache"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Let us perform intersection of reference peaks with span ones."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "import tempfile\n",
    "import pandas as pd\n",
    "from pandas import DataFrame\n",
    "\n",
    "def bed_length(cell_file):\n",
    "    t = pd.read_csv(cell_file,  sep='\\t', header=None)\n",
    "    if len(t) > 0:\n",
    "        return (t[2] - t[1]).sum()\n",
    "    return 0\n",
    "\n",
    "def lines(cell_file):\n",
    "    t = pd.read_csv(cell_file,  sep='\\t', header=None)\n",
    "    return len(t)\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",
    "def sorted_file(file):\n",
    "    ts = tempfile.mktemp()\n",
    "    ! cat {file} | sort -k1,1 -k2,2n -k3,3n > {ts}\n",
    "    return ts\n",
    "\n",
    "t = tempfile.mktemp()\n",
    "\n",
    "results = DataFrame(columns=['cell', 'fdr', 'ref_peaks', 'ref_length',\n",
    "                             'span_peaks', 'span_length', 'intersection_peaks', 'intersection_length'])\n",
    "for cell in tqdm(CELLS):\n",
    "    cell_file = sorted_file(f'{PATH}/peaks/{cell}_peaks_chr22.bed')\n",
    "    cell_peaks = lines(cell_file)\n",
    "    cell_peaks_length = bed_length(cell_file)    \n",
    "    \n",
    "    for reads in [f for f in listdir(PATH) if f.endswith(\".bam.bed\")]:\n",
    "        for fdr in FDRS: \n",
    "            file = f'{PATH}/span/{reads}_{fdr}_span.bed'\n",
    "            if cell not in file:\n",
    "                continue\n",
    "            ! cat {file} | grep chr22 > {t}\n",
    "            t = sorted_file(t)\n",
    "            span_peaks = lines(t)\n",
    "            span_length = bed_length(t)\n",
    "            intersection_file = f'{PATH}/peaks/intersection_{cell}_{fdr}.bed'\n",
    "            ! bedtools intersect -a {cell_file} -b {t} -wo > {intersection_file}\n",
    "            intersection_peaks = lines(intersection_file)\n",
    "            intersection_length = int(last_col(intersection_file).sum())\n",
    "            results.loc[len(results)] = [cell, fdr, cell_peaks, cell_peaks_length,\n",
    "                                         span_peaks, span_length, intersection_peaks, intersection_length]\n",
    "results"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have both reference peaks, span output and their intersection length, we can compute `TP`, `FP` and `FN` values."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "results['TP'] = results['intersection_length']\n",
    "results['FP'] = results['span_length'] - results['intersection_length']\n",
    "results['FN'] = results['ref_length'] - results['intersection_length']\n",
    "results['Spec'] = results['TP'] / (results['TP'] + results['FP'])\n",
    "results['Sen'] = results['TP'] / (results['TP'] + results['FN'])\n",
    "results.head(1)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "import seaborn as sns\n",
    "plt.figure(figsize=(4.5, 3))\n",
    "g_results = sns.lineplot(data=results,\n",
    "                         x='fdr', y='span_peaks', hue='cell',\n",
    "                         errorbar='se', marker='o', markersize=8)\n",
    "g_results.axes.title.set_text(f'Peaks')\n",
    "# g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "# plt.savefig(f'{GSE26320_PATH_HG38}/pics/H3K4me1_jaccard2.pdf', bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(4, 4))\n",
    "ax = plt.axes()\n",
    "for cell in tqdm(CELLS):\n",
    "    t = results[results['cell'] == cell].copy()\n",
    "    t.sort_values(by=['Sen', 'Spec'], inplace=True)\n",
    "    plt.plot(t['Sen'], t['Spec'], marker='o', label=cell)\n",
    "ax.set(xscale='log')\n",
    "plt.xlabel('Sensitivity')\n",
    "plt.ylabel('Specificity')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "CELLS = ['GM12878', 'HelaS3', 'K562']\n",
    "\n",
    "def Fb(b, spec, sen):\n",
    "    return (1 + b*b) * spec * sen / (b * b * spec + sen)\n",
    "\n",
    "plt.figure(figsize=(3, 4))\n",
    "ax = plt.axes()\n",
    "Beta = 0.5\n",
    "colors = {'GM12878': 'red', 'K562': 'blue', 'HelaS3':'violet'}\n",
    "plt.title('SPAN')\n",
    "plt.xlabel('threshold')\n",
    "plt.ylabel('F_0.5')\n",
    "for cell in CELLS:\n",
    "    r = results[(results['cell'] == cell)]\n",
    "    plt.plot(r['fdr'], Fb(Beta, r['Spec'], r['Sen']), linewidth=3.0, c=colors[cell], label=cell)\n",
    "plt.axvline(x=0.05, linewidth=3.0, c='grey', linestyle='--')\n",
    "plt.ylim(0.25, 0.85)\n",
    "# ax.set(xscale='log')\n",
    "plt.tight_layout()\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now we are ready to compare performance score F\n",
    "<img src=\"SPAN Comparison Opendnase.png\">"
   ]
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "language": "python",
   "display_name": "Python 3 (ipykernel)"
  },
  "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"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "source": [],
    "metadata": {
     "collapsed": false
    }
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}