Enh_Tlx_peaks_Mutation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "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",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load enhancers "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "enh = pb.BedTool('tracks/Enhancers_ChromHMM.bed')\n",
    "enh_df = enh.to_dataframe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load peaks and filter by score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tlx_pk = pb.BedTool('tracks/TLX3_TLX3_peaks.bed')\n",
    "tlx_df =tlx_pk.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": [
    "#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
}