build_nonoverlapping_TIFs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering non-overlapping transcripts for metagene analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pybedtools \n",
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### set script parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "home_dir = '/home/salma/Coding/cluster' \n",
    "\n",
    "#interval: genomic features are sepaced by at least 'interval' base pairs \n",
    "interval = 700\n",
    "clear_both_strands = False\n",
    "\n",
    "#input files\n",
    "TIF_input_file = os.path.join(home_dir, 'Annotations/steinmetz_transcripts_tifs_clean_sorted.gff')\n",
    "\n",
    "#output files\n",
    "TSS_nonoverlapping_file = os.path.join(home_dir, 'Annotations/salma/steinmetz_nonoverlapping_sense_TSS_' + str(interval) +'.gff')\n",
    "pA_nonoverlapping_file = os.path.join(home_dir, 'Annotations/salma/steinmetz_nonoverlapping_sense_pA_' + str(interval) +'.gff')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "tif = pd.read_table(TIF_input_file, sep='\\t', header=None)\n",
    "tif.columns = ['chr','ref','type','start','end','c1','strand','c2','gene_name']\n",
    "\n",
    "expanded_tif = tif.copy()\n",
    "expanded_tif.start = expanded_tif.start - int(interval/2)\n",
    "expanded_tif.end = expanded_tif.end + int(interval/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bed_expanded = pybedtools.bedtool.BedTool.from_dataframe(expanded_tif)\n",
    "\n",
    "def overlap_filter(feature):\n",
    "    return int(feature[3]) > 1\n",
    "\n",
    "bed_expanded_overlaps = bed_expanded.merge(c = 1, o = 'count').filter(overlap_filter)\n",
    "\n",
    "non_overlapping = bed_expanded.subtract(bed_expanded_overlaps, A = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "169"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#non-overlapping on both TSS and pA\n",
    "non_overlapping.count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#too little left\n",
    "#another strategy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### filtering overlaps on TSS and pA individually (here on sense strand)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "tif_bed = pybedtools.bedtool.BedTool.from_dataframe(tif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "upstream_TSS = list()\n",
    "for i in range(len(tif)):\n",
    "    if (tif.loc[i,'strand'] == '+'):\n",
    "        upstream_TSS.append([tif.loc[i,'chr'],\n",
    "                         '.','.',\n",
    "                         tif.loc[i,'start']-interval,\n",
    "                         tif.loc[i,'start']-1,\n",
    "                         i,'+','.', '.'])\n",
    "    else:\n",
    "        upstream_TSS.append([tif.loc[i,'chr'],\n",
    "                         '.','.',\n",
    "                         tif.loc[i,'end'] + 1,\n",
    "                         tif.loc[i,'end'] + interval,\n",
    "                         i,'-','.', '.'])\n",
    "upstream_TSS = pd.DataFrame(upstream_TSS, columns=tif.columns)\n",
    "upstream_TSS_bed = pybedtools.bedtool.BedTool.from_dataframe(upstream_TSS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3193"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TSS_nonoverlapping = upstream_TSS_bed.subtract(tif_bed, A = True, s = not clear_both_strands)\n",
    "\n",
    "indices = TSS_nonoverlapping.to_dataframe().loc[:,'score']\n",
    "TSS_nonoverlapping = tif.iloc[indices,:]\n",
    "pybedtools.bedtool.BedTool.from_dataframe(TSS_nonoverlapping).saveas(TSS_nonoverlapping_file)\n",
    "len(TSS_nonoverlapping)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "downstream_pA = list()\n",
    "for i in range(len(tif)):\n",
    "    if (tif.loc[i,'strand'] == '+'):\n",
    "        downstream_pA.append([tif.loc[i,'chr'],\n",
    "                         '.','.',\n",
    "                         tif.loc[i,'end'] + 1,\n",
    "                         tif.loc[i,'end'] + interval,\n",
    "                         i,'+','.', '.'])\n",
    "    else:\n",
    "        downstream_pA.append([tif.loc[i,'chr'],\n",
    "                         '.','.',\n",
    "                         tif.loc[i,'start'] - interval,\n",
    "                         tif.loc[i,'start'] - 1,\n",
    "                         i,'-','.', '.'])\n",
    "downstream_pA = pd.DataFrame(downstream_pA, columns=tif.columns)\n",
    "downstream_pA_bed = pybedtools.bedtool.BedTool.from_dataframe(downstream_pA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3193"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pA_nonoverlapping = downstream_pA_bed.subtract(tif_bed, A = True, s = not clear_both_strands)\n",
    "\n",
    "indices = pA_nonoverlapping.to_dataframe().loc[:,'score']\n",
    "pA_nonoverlapping = tif.iloc[indices,:]\n",
    "pybedtools.bedtool.BedTool.from_dataframe(pA_nonoverlapping).saveas(pA_nonoverlapping_file)\n",
    "len(pA_nonoverlapping)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}