SupplementalFigures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "39608876",
   "metadata": {},
   "source": [
    "# SFig. 3: HMM validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7db205f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def AIC(L, n_states):\n",
    "    n_params = 2 # because normal distribution defined by mean and standard deviation\n",
    "    return(2*(n_states**2 + n_params*n_states - 1) - 2*(L))\n",
    "\n",
    "def BIC(L, n_states, n_obs):\n",
    "    n_params = 2 # because normal distribution defined by mean and standard deviation\n",
    "    return((n_params*n_states + n_states**2 - 1)*np.log(n_obs) - 2*L)\n",
    "\n",
    "bics = []\n",
    "n_states_list = []\n",
    "aics = []\n",
    "cts = []\n",
    "\n",
    "for ix, row in org_df.query('chip == \"LB1\"').iterrows():\n",
    "    for n_states in [2, 3, 4, 5]:\n",
    "        ct = row['cell_type']\n",
    "        with open(f'/pollard/home/kathleen/projects/LADs/LADs_KDDs/hmms_no_umap_filter_{n_states}states/{ct}.json') as json_file:\n",
    "                    dat_model = HiddenMarkovModel.from_json(json.load(json_file))\n",
    "        \n",
    "        dat = pd.read_table(f'/pollard/home/kathleen/projects/LADs/LADs_KDDs/hmm_calls_no_umap_filter_{n_states}states/{ct}.tsv', names=['chrom','start','stop','score0','score1','hmm_pred'],\n",
    "                           header=0)\n",
    "\n",
    "        bic = BIC(dat_model.log_probability(dat[['score0','score1']].to_numpy()), n_states, len(dat))\n",
    "\n",
    "        bics.append(bic)\n",
    "\n",
    "        aic = AIC(dat_model.log_probability(dat[['score0','score1']].to_numpy()), n_states)\n",
    "\n",
    "        aics.append(aic)\n",
    "\n",
    "        n_states_list.append(n_states)\n",
    "\n",
    "        cts.append(ct)\n",
    "\n",
    "tog = pd.DataFrame({\n",
    "    'ct':cts,\n",
    "    'aic':aics,\n",
    "    'bic':bics,\n",
    "    'n_states':n_states_list\n",
    "}).drop_duplicates()\n",
    "\n",
    "aic_diffs = []\n",
    "bic_diffs = []\n",
    "states = []\n",
    "cts = []\n",
    "\n",
    "for state in [2,3,4]:\n",
    "    for ct in tog['ct'].unique():\n",
    "        cts.append(ct)\n",
    "        s_plus_one = state + 1\n",
    "        states.append(f'{state} > {s_plus_one}')\n",
    "        aic_diffs.append(tog.query('ct == @ct and n_states == @state')['aic'].tolist()[0]-tog.query('ct == @ct and n_states == @s_plus_one')['aic'].tolist()[0])\n",
    "        bic_diffs.append(tog.query('ct == @ct and n_states == @state')['bic'].tolist()[0]-tog.query('ct == @ct and n_states == @s_plus_one')['bic'].tolist()[0])\n",
    "        \n",
    "        \n",
    "\n",
    "df = pd.DataFrame({\n",
    "    'AIC_diff':aic_diffs,\n",
    "    'BIC_diff':bic_diffs,\n",
    "    'states':states,\n",
    "    'ct':cts\n",
    "})\n",
    "\n",
    "sns.barplot(x='ct', y='AIC_diff', hue='states', data=df,\n",
    "           order=sorted(df['ct'].drop_duplicates()))\n",
    "plt.xticks(rotation=45, ha='right')\n",
    "plt.title('LAD AIC differences: Models trained per cell type and per replicate')\n",
    "plt.legend(title='# of states', ncol=1, frameon=False)\n",
    "plt.ylabel('AIC difference')\n",
    "plt.xlabel('cell type')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "491a480d",
   "metadata": {},
   "source": [
    "# SFig. 5: KDD HMM validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b39eb1d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "bics = []\n",
    "n_states_list = []\n",
    "aics = []\n",
    "cts = []\n",
    "\n",
    "for ix, row in org_df.query('chip == \"H3K9me2\"').iterrows():\n",
    "    for n_states in [2, 3, 4, 5]:\n",
    "        ct = row['cell_type']\n",
    "        with open(f'/pollard/home/kathleen/projects/LADs/LADs_KDDs/hmms_no_umap_filter_{n_states}states/{ct}.json') as json_file:\n",
    "                    dat_model = HiddenMarkovModel.from_json(json.load(json_file))\n",
    "        \n",
    "        dat = pd.read_table(f'/pollard/home/kathleen/projects/LADs/LADs_KDDs/hmm_calls_no_umap_filter_{n_states}states/{ct}.tsv', names=['chrom','start','stop','score0','score1','hmm_pred'],\n",
    "                           header=0)\n",
    "\n",
    "        bic = BIC(dat_model.log_probability(dat[['score0','score1']].to_numpy()), n_states, len(dat))\n",
    "\n",
    "        bics.append(bic)\n",
    "\n",
    "        aic = AIC(dat_model.log_probability(dat[['score0','score1']].to_numpy()), n_states)\n",
    "\n",
    "        aics.append(aic)\n",
    "\n",
    "        n_states_list.append(n_states)\n",
    "\n",
    "        cts.append(ct)\n",
    "\n",
    "tog_kdd = pd.DataFrame({\n",
    "    'ct':cts,\n",
    "    'aic':aics,\n",
    "    'bic':bics,\n",
    "    'n_states':n_states_list\n",
    "}).drop_duplicates()\n",
    "\n",
    "aic_diffs = []\n",
    "bic_diffs = []\n",
    "states = []\n",
    "cts = []\n",
    "\n",
    "for state in [2,3,4]:\n",
    "    for ct in tog['ct'].unique():\n",
    "        cts.append(ct)\n",
    "        s_plus_one = state + 1\n",
    "        states.append(f'{state} > {s_plus_one}')\n",
    "        aic_diffs.append(tog.query('ct == @ct and n_states == @state')['aic'].tolist()[0]-tog.query('ct == @ct and n_states == @s_plus_one')['aic'].tolist()[0])\n",
    "        bic_diffs.append(tog.query('ct == @ct and n_states == @state')['bic'].tolist()[0]-tog.query('ct == @ct and n_states == @s_plus_one')['bic'].tolist()[0])\n",
    "        \n",
    "        \n",
    "\n",
    "kdd_df = pd.DataFrame({\n",
    "    'AIC_diff':aic_diffs,\n",
    "    'BIC_diff':bic_diffs,\n",
    "    'states':states,\n",
    "    'ct':cts\n",
    "})\n",
    "\n",
    "ct_replace = {\n",
    "    'CardiacMyocytes':'Cardiac myocytes',\n",
    "    'EndoProgenitor':'Endothelial progenitors',\n",
    "    'D5Midbrain':'Midbrain progenitors',\n",
    "    'H9ESC':'ESCs',\n",
    "    'MidHindgut':'Mid-hindgut',\n",
    "    'ParaxMesoderm':'Paraxial mesoderm',\n",
    "    'Early Somite':'Early somite',\n",
    "    'DefEctoderm':'Definitive ectoderm',\n",
    "    'BorderEctoderm':'Border ectoderm',\n",
    "    'D4Artery':'Artery progenitors',\n",
    "    'APS':'Anterior primitive streak',\n",
    "    'EarlySomite':'Early somite'\n",
    "}\n",
    "\n",
    "sns.barplot(x='ct', y='AIC_diff', hue='states', data=kdd_df,\n",
    "           order=sorted(df['ct'].drop_duplicates()))\n",
    "plt.xticks(rotation=45, ha='right')\n",
    "plt.title('KDD AIC differences: Models trained per cell type and per replicate')\n",
    "plt.legend(title='# of states', ncol=1, frameon=False)\n",
    "plt.ylabel('AIC difference')\n",
    "plt.xlabel('cell type')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de23c00c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# H3K9me2 in KDDs\n",
    "\n",
    "dat = pd.read_csv('tables_out/all_dat_tog.tsv', sep='\\t')\n",
    "\n",
    "dat['score_lb1'] = dat[['score0_lb1','score1_lb1']].mean(axis=1)\n",
    "dat['score_k9'] = dat[['score0_k9','score1_k9']].mean(axis=1)\n",
    "\n",
    "c = sns.violinplot(y='score_k9', x='cell_type', data=dat, kind='violin',\n",
    "           hue='category_k9', hue_order=['nonKDD','T2-KDD','T1-KDD'],\n",
    "                   order=sorted(dat['cell_type'].unique()),\n",
    "              linewidth=0.5, inner=None, ax=axC)\n",
    "plt.xticks(rotation=45, ha='right')\n",
    "plt.title('H3K9me2 binding, KDDs')\n",
    "plt.ylabel('H3K9me2 binding [log2(read counts)]')\n",
    "plt.xlabel('cell type')\n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "axC.legend(ncol=3, title='LAD category', frameon=False, bbox_to_anchor=(0.98, 1.05))\n",
    "axC.set(ylim=(-3, 3.8))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "617ddddf",
   "metadata": {},
   "source": [
    "# SFig. 6: KDD/LAD overlap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "652f054a",
   "metadata": {},
   "outputs": [],
   "source": [
    "cts = ['Liver', 'CardiacMyocytes', 'EndoProgenitor', 'D5Midbrain',\n",
    "       'H9ESC', 'Epicardium', 'MidHindgut', 'ParaxMesoderm',\n",
    "       'EarlySomite', 'DefEctoderm', 'BorderEctoderm', 'D4Artery', 'APS']\n",
    "\n",
    "t1lads = {}\n",
    "t1kdds = {}\n",
    "t2lads = {}\n",
    "t2kdds = {}\n",
    "\n",
    "for ct in cts:\n",
    "    t1lads[ct] = BedTool(f'{path_to_beds}/BED_files_no_umap_filter_3states/{ct}_T1-LADs.bed').sort()\n",
    "    t1kdds[ct] = BedTool(f'{path_to_beds}/KDD_BED_files_no_umap_filter_3states/{ct}_T1-KDDs.bed').sort()\n",
    "    t2lads[ct] = BedTool(f'{path_to_beds}/BED_files_no_umap_filter_3states/{ct}_T2-LADs.bed').sort()\n",
    "    t2kdds[ct] = BedTool(f'{path_to_beds}/KDD_BED_files_no_umap_filter_3states/{ct}_T2-KDDs.bed').sort()\n",
    "    \n",
    "\n",
    "f, ax = plt.subplots(figsize=(4.5,3))\n",
    "\n",
    "lad_types = []\n",
    "percs = []\n",
    "\n",
    "for ct in cts:\n",
    "    t1kdd = t1kdds[ct]\n",
    "    t1lad = t1lads[ct]\n",
    "    t2kdd = t2kdds[ct]\n",
    "    t2lad = t2lads[ct]\n",
    "    percs.append(100.0*len(t1lad.intersect(t1kdd, u=True))/len(t1lad))\n",
    "    lad_types.append('T1-LAD/KDD')\n",
    "    percs.append(100.0*len(t2lad.intersect(t2kdd, u=True))/len(t2lad))\n",
    "    lad_types.append('T2-LAD/KDD')\n",
    "\n",
    "sim_plot = pd.DataFrame({\n",
    "    'percentage':percs,\n",
    "    'type':lad_types\n",
    "})\n",
    "\n",
    "sns.barplot(x='type', y='percentage', data=sim_plot, ci='sd',\n",
    "           palette=dict(zip(sim_plot['type'].unique(),['darkgrey','dimgrey'])))\n",
    "plt.ylim(0, 100)\n",
    "plt.ylabel('LADs overlapping KDDs (%)')\n",
    "plt.xlabel('Type')\n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\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
}