Enh_Genes_Tlx_peaks.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pybedtools as pb\n",
    "import numpy as np\n",
    "#import scipy\n",
    "import pandas as pd\n",
    "#import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "%load_ext autoreload \n",
    "%autoreload 2  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from os.path import join \n",
    "WORKDIR = '/home/sergio/Res_CIML/TLX3_project'\n",
    "SCRIPTS = join(WORKDIR,'scripts')\n",
    "DATADIR = join(WORKDIR,'data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load enhancers "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "enh = pb.BedTool(join(DATADIR,'tracks/Enhancers_ChromHMM.bed'))\n",
    "enh_df = enh.to_dataframe()\n",
    "\n",
    "# Enh to genes table\n",
    "enh2gn = pd.read_table(join(DATADIR,'tracks/Enhancers_ChromHMM_enh2genes.txt'), \n",
    "                        header=1, \n",
    "                        names=['enhancers','gene_name'])\n",
    "# Gene to enh table\n",
    "gn2enh = pd.read_table(join(DATADIR,'tracks/Enhancers_ChromHMM_genes2enh.txt'), \n",
    "                        header=1, \n",
    "                        names=['gene_name','enhancers'])\n",
    "\n",
    "def enh_gene(genl,df):\n",
    "    if isinstance(genl, str):\n",
    "        genl = [genl]\n",
    "    enhll = list()\n",
    "    for gen in genl:\n",
    "        reg = df[df['gene_name']==gen]\n",
    "        if len(reg)>0:\n",
    "            dr = reg.iloc[0]['enhancers'].split(' ')\n",
    "            enhl = [x for x in dr if 'enh' in x]\n",
    "        else:\n",
    "            enhl=[]\n",
    "        enhll = enhll+enhl\n",
    "    return list(set(enhll))\n",
    "\n",
    "def gene_enh(enhl,df):\n",
    "    if isinstance(enhl, str):\n",
    "        enhl = [enhl]\n",
    "    \n",
    "    genl = list()\n",
    "    for enh in enhl:\n",
    "        reg = df[df['enhancers']==enh]\n",
    "        if len(reg)>0:\n",
    "            ls = reg.iloc[0]['gene_name'].split(', ')\n",
    "            gl=[tr.split(' (')[0] for tr in ls]\n",
    "        else:\n",
    "            gl=[]\n",
    "        genl = genl+gl\n",
    "    return list(set(genl))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['enh_6659', 'enh_18785', 'enh_6657', 'enh_18786', 'enh_18790', 'enh_18777', 'enh_6658', 'enh_18779', 'enh_18781', 'enh_18787', 'enh_18783', 'enh_6661', 'enh_6660', 'enh_18791', 'enh_18789', 'enh_18793', 'enh_18778', 'enh_18780', 'enh_18792', 'enh_18794', 'enh_18784', 'enh_18788', 'enh_6662', 'enh_18782']\n",
      "['Il1r2', 'Map4k4']\n"
     ]
    }
   ],
   "source": [
    "gens = ['Tm9sf2', 'Tle3']\n",
    "gnn = 'Tm9sf2'\n",
    "enhs = ['enh_168', 'enh_1683', 'enh_268']\n",
    "enhh = 'enh_168'\n",
    "\n",
    "print(enh_gene(gens,gn2enh))\n",
    "print(gene_enh(enhh,enh2gn))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load peaks "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tlx = pb.BedTool('tracks/TLX3_TLX3_peaks.bed')\n",
    "tlx_df =tlx.to_dataframe(names=['chr','start','end','name','score'])\n",
    "\n",
    "#shr = 100\n",
    "#top = tlx_df[tlx_df['score']>shr].sort_values('score', axis=0, ascending=False)\n",
    "#tlx_top = pb.BedTool.from_dataframe(top)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#(enh+tlx).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load TSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "genes = pb.BedTool('tracks/annot_tracks/genes.bed')\n",
    "\n",
    "### == TSS region\n",
    "tss =  3000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gen_wth_peaks(en,enh2gn,tlx,genes,tss=3000):\n",
    "    gn_tss = genes.slop(b=tss, genome='mm9')\n",
    "    gn_tss_tlx=gn_tss+tlx\n",
    "    gn_tss_tlx_df = gn_tss_tlx.to_dataframe()\n",
    "    \n",
    "    gn_tlx_l =list(gn_tss_tlx_df['name']) \n",
    "    \n",
    "    gl = gene_enh(en,enh2gn)\n",
    "    \n",
    "    return list(set(gl) & set(gn_tlx_l))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "enh_tlx = enh+tlx\n",
    "enh_tlx_df = enh_tlx.to_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "enh_tlx_df['gentlx'] = enh_tlx_df.apply(lambda row: gen_wth_peaks(row['name'],enh2gn,tlx,genes), axis=1)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "enh_tlx_df['num_gentlx'] = enh_tlx_df.apply(lambda row: len(row['gentlx']), axis=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "enh_tlx_df_sort = enh_tlx_df[enh_tlx_df['num_gentlx']>0]\n",
    "enh_tlx_df_sort.sort_values('num_gentlx', axis=0, ascending=False,inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#enh_tlx_df_sort.head()\n",
    "len(enh_tlx_df_sort)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load mutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ftlx = 'tracks/WGS-WES/Germline/FERRIER_09_Germline.allchr.snpEff.p.SAL.SAL10_1.vcf'\n",
    "var_b = pb.BedTool(ftlx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#allTLX_list =list(all_tss_gene['name'])\n",
    "#allTLX_list.sort()\n",
    "#len(allTLX_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#a,b,c = plt.hist(np.log(tlx_df['score']), bins=100)\n",
    "#a,b,c = plt.hist(tlx_df['score'], bins=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intersect enhancers and high score peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#enh_tlx_top = enh+tlx_top "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#enh_tlx_top.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load WGS variants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "path = 'tracks/WGS-WES/Germline/'\n",
    "ftlx = 'tracks/WGS-WES/Germline/FERRIER_09_Germline.allchr.snpEff.p.SAL.SAL10_1.vcf'\n",
    "var_b = pb.BedTool(ftlx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find variants in enhancers with TLX3 peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "var_ehn_tlx = var_b+enh_tlx_top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "var_ehn_tlx_df = var_ehn_tlx.to_dataframe(names=['CHROM','POS','ID','REF','ALT','QUAL','FILTER','INFO','FORMAT','AC3812','AC3813','AC3814','AC3815'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#var_ehn_tlx_df.head(10)\n",
    "len(var_ehn_tlx_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Enhancers Table manipulation\n",
    "enh2gn = pd.read_table('tracks/Enhancers_ChromHMM_enh2genes.txt', \n",
    "                        header=1, \n",
    "                        names=['enhancers','gene_name'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "enh_tlx_top_df = enh_tlx_top.to_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "enh_nm = list(enh_tlx_top_df['name'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "enh2gn_top = enh2gn[enh2gn['enhancers'].isin(enh_nm)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = list()\n",
    "for i in range(len(enh2gn_top )):\n",
    "    ls = enh2gn_top['gene_name'].iloc[i].split(', ')\n",
    "    a=a+[tr.split(' (')[0] for tr in ls]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "len(set(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Enrichment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "import EnrichRLib as erl\n",
    "# List of gene sets as above\n",
    "gss = [ \n",
    "       'GO_Biological_Process_2018',\n",
    "       'GO_Cellular_Component_2018',\n",
    "       'GO_Molecular_Function_2018',\n",
    "       'KEGG_2016',\n",
    "       'Reactome_2016'\n",
    "       ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = [x.upper() for x in a]\n",
    "\n",
    "enr_a = erl.enrich_gs(a,gss)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# --- Plot ---\n",
    "enr_a.sort_values('p-Val', axis=0, inplace = True)\n",
    "ds = enr_a.head(20)\n",
    "\n",
    "f, ax = plt.subplots()\n",
    "sns.barplot(y=ds.index,\n",
    "            x='-log10(p-Val)',\n",
    "            ax = ax, \n",
    "            color=\"Red\", \n",
    "            data = ds)\n",
    "ax.set_title('Genes regulated by enhancers with high Tlx3 peaks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Expression analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import RNA_expression_processing as rn\n",
    "from os.path import join"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tbl = pd.read_table(join('tracks', 'TLX3vsRAG-results_genesNames.txt'), index_col=0)\n",
    "#tbl = tbl[(tbl.padj < 0.05)].dropna()\n",
    "\n",
    "# === Pheno ==\n",
    "A,B = 'TLX3','RAG'\n",
    "classes = [A]*3+[B]*3\n",
    "\n",
    "\n",
    "cols = ['Gene_name', 'TLX3.1_1','TLX3.1_5','TLX3.1_P','R2.RAG1W.RAG1','RAGS.RAGZ','RAGZ']\n",
    "\n",
    "tbn = tbl[cols]\n",
    "tbv = tbn.set_index(keys=tbn.columns[0])\n",
    "tbv.index=tbv.index.str.upper()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "### == UP analysis\n",
    "\n",
    "tbu = tbv.iloc[tbv.index.isin(a)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "jj,kk,ll=rn.scatter_n(tbu, A, B, classes, n_top=3, ttl='') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ntop=40\n",
    "gr = rn.cluster(tbu, A, B, classes, n_top=ntop)\n",
    "gr.ax_heatmap.set_title('Cluster '+A+'/'+B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.5.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}