Figure1_panels.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "207035de",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib as mpl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2d14300",
   "metadata": {},
   "source": [
    "# LB1 density plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2da47787",
   "metadata": {},
   "outputs": [],
   "source": [
    "density_dat = pd.read_csv('tables_out/all_dat_tog.tsv', sep='\\t').replace('EarlySomite','Early somite')\n",
    "density_dat['score_lb1'] = density_dat[['score0_lb1','score1_lb1']].mean(axis=1)\n",
    "density_dat['score_k9'] = density_dat[['score0_k9','score1_k9']].mean(axis=1)\n",
    "\n",
    "mypal = {'T1-LAD':'indigo','T2-LAD':'mediumpurple','nonLAD':'lightgrey'}\n",
    "\n",
    "f, axB = plt.subplots(figsize=(4,8))\n",
    "\n",
    "b = sns.violinplot(x='score_lb1', y='cell_type', data=density_dat, kind='violin',\n",
    "           hue='category_lb1',linewidth=0.5, inner=None, ax=axB,\n",
    "                   hue_order=['T1-LAD','T2-LAD','nonLAD'],\n",
    "          order=['ESCs','Mid-hindgut','Liver','Anterior primitive streak',\n",
    "                'Paraxial mesoderm','Early somite','Artery progenitors',\n",
    "                'Endothelial progenitors','Cardiac myocytes','Epicardium',\n",
    "                'Definitive ectoderm','Border ectoderm','Midbrain progenitors'],\n",
    "                  palette = mypal)\n",
    "plt.legend(frameon=False, title='LAD category', bbox_to_anchor=(1, 1.15),\n",
    "          ncol=1)\n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "axB.set(ylabel='cell type', xlabel='LB1 binding [log2(read counts)]')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "800bd321",
   "metadata": {},
   "source": [
    "# LB1 at LAD boundaries\n",
    "\n",
    "Example deepTools commands:\n",
    "\n",
    "`computeMatrix reference-point -S bigwigs/H9ESC_rep1.bigwig -R ctcf_analysis/non_to_t2_escs.bed --referencePoint center -a 500000 -b 500000 -out LB1_non_to_t2_escs.tab.gz -bs 20000\n",
    "\n",
    "\n",
    "plotProfile -m LB1_non_to_t2_escs.tab.gz --outFileName LB1_non_to_t2_escs_profile.png --outFileNameData LB1_non_to_t2_escs_for_plot.tsv`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6538738",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_deeptools_dat(path, rtype, dtype):\n",
    "    dat = pd.read_table(path,\n",
    "                         skiprows=1, header=0, low_memory=False).replace(np.nan, 0.0).replace('--',0.0)\n",
    "    dat.columns = [dtype ,rtype] + list(range(dat.shape[1]))[:-2]\n",
    "    dat.columns = [dtype,rtype] + list(range(dat.shape[1]))[:-2]\n",
    "\n",
    "    dfs = []\n",
    "\n",
    "    for ix, row in dat.iterrows():\n",
    "        dset = row[dtype]\n",
    "        r_dset = row[rtype]\n",
    "\n",
    "        dat = pd.DataFrame(row[2:].T)\n",
    "        dat.columns=[f'{dtype}_sig']\n",
    "        dat[dtype] = dset.split('_')[0]\n",
    "        dat[rtype] = r_dset.split('_')[0]\n",
    "        dat['bin'] = dat.index\n",
    "\n",
    "        dfs.append(dat)\n",
    "\n",
    "\n",
    "    tog = pd.concat(dfs, sort=False)\n",
    "    tog[f'{dtype}_sig'] = tog[f'{dtype}_sig'].astype(float)\n",
    "    return(tog)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0379c804",
   "metadata": {},
   "outputs": [],
   "source": [
    "# deepTools has a weird quirk where it outputs additional bins\n",
    "\n",
    "max_bin = 50\n",
    "\n",
    "# T2 to non\n",
    "\n",
    "non_to_t2_all = []\n",
    "\n",
    "for ct in ['escs','parax_mesoderm','cms']:\n",
    "    t2_to_non = load_deeptools_dat(f'LB1_t2_to_non_{ct}_for_plot.tsv','T2_to_non','LB1')#.query('bin > 0')\n",
    "    rev_bins = t2_to_non['bin'].tolist()\n",
    "    rev_bins.reverse()\n",
    "    t2_to_non['bin'] = rev_bins\n",
    "    non_to_t2 = load_deeptools_dat(f'LB1_non_to_t2_{ct}_for_plot.tsv','T2_to_non','LB1')#.query('bin < @max_bin')\n",
    "    t2_to_non_tog = pd.concat([t2_to_non, non_to_t2], sort=False).sort_values(by='bin', ascending=False)\n",
    "    t2_to_non_tog['cat'] = 'T2_to_non'\n",
    "    t2_to_non_tog['ct'] = ct\n",
    "    non_to_t2_all.append(t2_to_non_tog)\n",
    "    \n",
    "non_to_t2_all_cts = pd.concat(non_to_t2_all, sort=False)\n",
    "\n",
    "# T2 to T1\n",
    "\n",
    "t2_to_t1_all = []\n",
    "\n",
    "for ct in ['escs','parax_mesoderm','cms']:\n",
    "    t1_to_t2 = load_deeptools_dat(f'LB1_t1_to_t2_{ct}_for_plot.tsv','T2_to_T1','LB1')#.query('bin < @max_bin')\n",
    "    rev_bins = t1_to_t2['bin'].tolist()\n",
    "    rev_bins.reverse()\n",
    "    t1_to_t2['bin'] = rev_bins\n",
    "    t2_to_t1 = load_deeptools_dat(f'LB1_t2_to_t1_{ct}_for_plot.tsv','T2_to_T1','LB1')#.query('bin < @max_bin')\n",
    "    t2_to_t1_tog = pd.concat([t2_to_t1, t1_to_t2], sort=False).sort_values(by='bin', ascending=False)\n",
    "    t2_to_t1_tog['cat'] = 'T2_to_t1'\n",
    "    t2_to_t1_tog['ct'] = ct\n",
    "    t2_to_t1_all.append(t2_to_t1_tog)\n",
    "    \n",
    "t2_to_t1_all_cts = pd.concat(t2_to_t1_all, sort=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbba4c2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,3))\n",
    "gs = mpl.gridspec.GridSpec(1, 2, width_ratios=[1, 1],\n",
    "                          height_ratios=[1])\n",
    "\n",
    "# nonLAD --> T2-LAD\n",
    "\n",
    "axA = plt.subplot(gs[0, 0])\n",
    "a = sns.lineplot(x='bin', y='LB1_sig', data=non_to_t2_all_cts.replace({'escs':'ESCs',\n",
    "                                                                       'cms':'Cardiac myocytes',\n",
    "                                                                       'parax_mesoderm':'Paraxial mesoderm'}), \n",
    "                 hue='ct', ci='sd', alpha=0.5, estimator='mean', \n",
    "                 hue_order=['ESCs','Paraxial mesoderm','Cardiac myocytes'],\n",
    "                ax=axA)\n",
    "a.set(ylabel='LB1 [log2(read counts)]', xlabel='genomic position', ylim=(-0.4, 0.7))\n",
    "axA.legend(frameon=False, bbox_to_anchor = (1., -0.2))\n",
    "axA.set_xticks([0, 25, 50])\n",
    "axA.set_xticklabels(['-0.5Mb','boundary','+0.5Mb'])\n",
    "axA.axvline(25, linestyle='--', color='grey')\n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "axA.annotate('',\n",
    "            xy=(0.1, .98), xycoords='axes fraction',\n",
    "            xytext=(0.4, .98), textcoords='axes fraction',\n",
    "            arrowprops=dict(facecolor='grey', shrink=0.05, alpha=0.5),\n",
    "            horizontalalignment='right', verticalalignment='top')\n",
    "axA.annotate('',\n",
    "            xy=(0.9, .98), xycoords='axes fraction',\n",
    "            xytext=(0.6, .98), textcoords='axes fraction',\n",
    "            arrowprops=dict(facecolor='grey', shrink=0.05, alpha=0.5),\n",
    "            horizontalalignment='right', verticalalignment='top')\n",
    "\n",
    "# T2-LAD --> T1-LAD\n",
    "\n",
    "axB = plt.subplot(gs[0, 1])\n",
    "b = sns.lineplot(x='bin', y='LB1_sig', data=t2_to_t1_all_cts.replace({'escs':'ESCs',\n",
    "                                                                       'cms':'Cardiac myocytes',\n",
    "                                                                       'parax_mesoderm':'Paraxial mesoderm'}), \n",
    "                 hue='ct', ci='sd', alpha=0.5, estimator='mean', \n",
    "                 hue_order=['ESCs','Paraxial mesoderm','Cardiac myocytes'],\n",
    "                ax=axB)\n",
    "b.set(ylabel='LB1 [log2(read counts)]', xlabel='genomic position', ylim=(-0.4, 0.7))\n",
    "axB.legend(frameon=False, bbox_to_anchor = (1., -0.2))\n",
    "axB.set_xticks([0, 25, 50])\n",
    "axB.set_xticklabels(['-0.5Mb','boundary','+0.5Mb'])\n",
    "axB.axvline(25, linestyle='--', color='grey')\n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "axB.annotate('',\n",
    "            xy=(0.1, .98), xycoords='axes fraction',\n",
    "            xytext=(0.4, .98), textcoords='axes fraction',\n",
    "            arrowprops=dict(facecolor='grey', shrink=0.05, alpha=0.5),\n",
    "            horizontalalignment='right', verticalalignment='top')\n",
    "axB.annotate('',\n",
    "            xy=(0.9, .98), xycoords='axes fraction',\n",
    "            xytext=(0.6, .98), textcoords='axes fraction',\n",
    "            arrowprops=dict(facecolor='grey', shrink=0.05, alpha=0.5),\n",
    "            horizontalalignment='right', verticalalignment='top')\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LAD_analysis",
   "language": "python",
   "name": "lad_analysis"
  },
  "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}