hepatocytes.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sustained-diving",
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "from IPython.display import display\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import warnings\n",
    "import pickle\n",
    "%matplotlib inline\n",
    "#%matplotlib tk\n",
    "'''\n",
    "import logging\n",
    "\n",
    "logging.basicConfig(filename='c14.log',\n",
    "                              filemode='a',\n",
    "                              format='%(asctime)s,%(msecs)d %(name)s %(levelname)s %(message)s',\n",
    "                             datefmt='%H:%M:%S',\n",
    "                             level=logging.DEBUG)\n",
    "'''                 \n",
    "pd.set_option('display.max_rows', 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "competitive-bundle",
   "metadata": {},
   "outputs": [],
   "source": [
    "import arviz as az\n",
    "from matplotlib import rcParams,rc\n",
    "import seaborn as sb\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "fonts = 1\n",
    "sb.set(context=\"paper\",style='ticks',font_scale=fonts)\n",
    "plt.rcParams[\"font.family\"] = \"Times New Roman\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sublime-biography",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import c14\n",
    "import c14.models.liver_new_mitosis as lm\n",
    "import c14.corner\n",
    "from c14.utils import *\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "european-samuel",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "necessary-demographic",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Data files can be found in Zenodo \n",
    "directory = 'hepatocytes/'\n",
    "data = pd.read_csv('C14data_liver_samples_20211020.csv')\n",
    "data = data.groupby(['type', 'sample', 'ploidy', 'pathology']).mean().dropna(how='all').reset_index()\n",
    "data['age'] = data['Dcoll'] - data['Dbirth']\n",
    "data = data.query('type == \"hepatocyte\" and pathology != [\"Y\", \"C\"]')\n",
    "data['measurment_types'] = data['ploidy']\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "composite-gather",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tired-administrator",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.seterr(all='ignore')\n",
    "results= dict()\n",
    "NUM=0\n",
    "results= dict(data=data,results={},models={})\n",
    "\n",
    "\n",
    "for mm in lm.models_list :\n",
    "    m = mm(dnatotal=lm.dnatotal_spline)\n",
    "    if m.__class__.__name__ not in ['POP1','Rm','POP3p_2x2n_NP']:\n",
    "        continue\n",
    "        pass\n",
    "    path = directory+m.__class__.__name__+'spline_'+str(NUM) + '.pickle'\n",
    "    try:\n",
    "        with open(path, 'rb') as handle:\n",
    "            res = pickle.load(handle)\n",
    "    except FileNotFoundError:\n",
    "        print(m.__class__.__name__,path)\n",
    "        continue\n",
    "    except:\n",
    "        print(path)\n",
    "        continue\n",
    "    #print('FOUND',mm.__name__)\n",
    "    results['results'].update(  {m.__class__.__name__+'spline' : {'raw':res} } )\n",
    "    results['models'].update(   {m.__class__.__name__+'spline' : m } )        \n",
    "\n",
    "\n",
    "get_arviz(results,burnin=1000)\n",
    "ranking = get_getranking(results)\n",
    "ranking\n",
    "printt = get_pointestimates(results,error_in_real=False)[['Parameter','Value','Confidence Interval']]\n",
    "printt.to_excel('liver_fix_is_newploidy.xlsx')\n",
    "calc_c14_all(results)\n",
    "printt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "great-layout",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import normaltest\n",
    "from scipy.stats import shapiro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "personalized-identification",
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_val = {i :[] for i in results['models'].keys()}\n",
    "for i,(model_name,model) in enumerate(results['models'].items()):\n",
    "    results_m = results['results'][model_name]\n",
    "    for subject,row in data.iterrows():\n",
    "        sigf = results['point_est'].loc[model_name,('median','sigma')]\n",
    "        delta = (results_m['c14_M'].loc[subject]-row['d14C'])/np.sqrt(row['e14C']**2+sigf**2)\n",
    "        if np.abs(delta)<0.1 or True:\n",
    "            delta_val[model_name].append(delta)\n",
    "    print(model_name,normaltest(delta_val[model_name]))\n",
    "    print(model_name,shapiro(delta_val[model_name]))\n",
    "    \n",
    "    plt.hist(delta_val[model_name],density=True)\n",
    "    plt.title(model_name)\n",
    "    plt.figure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fancy-barcelona",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "\n",
    "subject = 1\n",
    "excel_writer = pd.ExcelWriter('time_rates_flows_newploidy.xlsx')\n",
    "for m_i,(model_name,model) in enumerate(results['models'].items()):\n",
    "    print(model_name)\n",
    "    age = results['data'].loc[subject,'age']\n",
    "    T = np.arange(0,90,0.1)\n",
    "    parameter_phy =  {i:results['point_est'].loc[model_name].loc[('median',i)] for i in model.parameter_names}\n",
    "    model.set_parameters_phy(parameter_phy,mode='bayes')\n",
    "    ipsd=[]\n",
    "    for i in T:\n",
    "        ip = model.calc_implicit_parameters(i)\n",
    "        ipsd.append(ip)\n",
    "\n",
    "    a_i =  [j for j,x in enumerate(results['data'].index) if x == subject][0]  \n",
    "    ipid = listofdict_to_dictofarray_f(ipsd,a_i)\n",
    "    ipid.update({n:np.array([v]*len(T)) for n,v in parameter_phy.items()})\n",
    "    df = pd.DataFrame(ipid)\n",
    "    df.index.name = 'age'\n",
    "    df.index =df.index/10\n",
    "    \n",
    "    \n",
    "    flow_name = defaultdict(list)\n",
    "    for i,(ind2,flows) in enumerate(model.flow_in.items()):\n",
    "        y = np.zeros_like(T)\n",
    "        for rate,pop,factor in flows:\n",
    "            try:\n",
    "                if pop is None:\n",
    "                    y += ipid[rate] * factor\n",
    "                else:\n",
    "                    y += ipid[rate] * factor * ipid[pop]\n",
    "            except Exception as e:\n",
    "                raise Exception([e,Exception(f'in class {model_name}')])\n",
    "        #df['inflow_'+ind2] = y\n",
    "        for rate,pop,factor in flows:\n",
    "            try:\n",
    "                if pop is None:\n",
    "                    yi = ipid[rate] * factor\n",
    "                else:\n",
    "                    yi = ipid[rate] * factor * ipid[pop]\n",
    "            except Exception as e:\n",
    "                raise Exception([e,Exception(f'in class {model_name}')])\n",
    "            df['probIN_'+pop+'_to_'+ind2] = yi/y\n",
    "            df['flow_'+pop+'_to_'+ind2] = yi\n",
    "            flow_name[pop] += [ind2]\n",
    "        #print(ind2,flows)\n",
    "    for pop_source,val in flow_name.items():\n",
    "        norm = 0\n",
    "        for pop in val:\n",
    "            norm += df['flow_'+pop_source+'_to_'+pop]\n",
    "        for pop in val:\n",
    "            df['probOUT_'+pop_source+'_to_'+pop] = df['flow_'+pop_source+'_to_'+pop]/norm\n",
    "    df.to_excel(excel_writer=excel_writer,sheet_name=model_name)\n",
    "    if model_name =='Rmspline':\n",
    "        print(df[['b2n','b4n']].loc[50]*np.array([model.ploidy(50),1-model.ploidy(50)])*1e11*1/365/1e6)\n",
    "        N50 = np.sum(df[['N4n','N2n']].loc[50])\n",
    "        print('N2',np.sum(df[['flow_N2n_to_N2n','flow_N4n_to_N2n']].loc[50])/N50*1e11*1/365/1e6)\n",
    "        print('N4',np.sum(df[['flow_N2n_to_N4n','flow_N4n_to_N4n']].loc[50])/N50*1e11*1/365/1e6)\n",
    "    if model_name =='POP3p_2x2n_NPspline':\n",
    "        print(df[['g2nb2n','k2x2nb2x2n','k4nb4n']].loc[50]*np.array([model.ploidy(50),model.ploidy2x2(50),1-model.ploidy2x2(50)-model.ploidy(50)])*1e11*1/365/1e6)\n",
    "        N50 = np.sum(df[['N4n','N2x2n','N2n']].loc[50])\n",
    "        NN50 = 2.5225*1e11\n",
    "        print('N2',np.sum(df[['flow_N2x2n_to_N2n','flow_N2n_to_N2n','flow_N4n_to_N2n']].loc[50])/N50*NN50*1/365/1e6)\n",
    "        print('N2x2',np.sum(df[['flow_N2x2n_to_N2x2n','flow_N4n_to_N2x2n']].loc[50])/N50*NN50*1/365/1e6)\n",
    "        print('N4',np.sum(df[['flow_N2x2n_to_N4n','flow_N4n_to_N4n']].loc[50])/N50*NN50*1/365/1e6)\n",
    "\n",
    "excel_writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "grand-madness",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "annoying-tuner",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "antique-casting",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "seven-vietnam",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "friendly-century",
   "metadata": {},
   "outputs": [],
   "source": [
    "with pd.HDFStore('point_est.pandas') as st:\n",
    "    st['F_IS_spline_newploidy'] = results['point_est']\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "corporate-vinyl",
   "metadata": {},
   "outputs": [],
   "source": [
    "rename={'POP3p_2x2n_NPspline':{'sigma':r'$\\sigma$','d2n':r'$\\delta_2$','d2x2n':r'$\\delta_{2x2}$','d4n':r'$\\detla_4$',\n",
    "                               'k2nb4n':r'$\\kappa_{4,2}$','g2x2nb2n':r'$\\kappa_{2x2,2}$','k2nb2x2n':r'$\\kappa_{2,2x2}$'},\n",
    "        'Rmspline':{'sigma':r'$\\sigma$','d2n':r'$\\delta_2$','d4n':r\"$\\delta_{p}$\",'k24':r\"$\\kappa_{2p}$\",'k42':r\"$\\kappa_{p2}$\"},\n",
    "        'POP1spline':{'sigma':r'$\\sigma$','d':r'$\\delta$'}}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "oriented-voluntary",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import matplotlib.patches as patches\n",
    "import matplotlib as mpl\n",
    "for n,v in results['results'].items():  \n",
    "    fimag = 'overviewF_IS_'+n+'.png'\n",
    "    if not os.path.isfile(fimag) or True :\n",
    "        c14.corner.corner_R(v['azdata'].posterior.a.values,\n",
    "                   results['models'][n].parameter_names,\n",
    "                        point_estimate=results['models'][n].transform_physical_to_fit(v['median'],mode='bayes'),\n",
    "                        rename=rename[n],unitlog='[1/year]',\n",
    "                        logparas=results['models'][n].logparas)    \n",
    "        f = plt.gcf()\n",
    "        f.show();\n",
    "        if n=='POP1spline':\n",
    "            bbax0 = f.get_axes()[0].get_position()\n",
    "            bbax3 = f.get_axes()[2].get_position()\n",
    "            ticks = [-1,0]\n",
    "            f.get_axes()[0].set_xticks(ticks)\n",
    "            f.get_axes()[0].set_xticklabels([r\"$10^{\"+str(int(i))+r\"}$\"  for i in ticks])\n",
    "            \n",
    "            f.get_axes()[1].set_xticks(ticks)\n",
    "            f.get_axes()[1].set_xticklabels([r\"$10^{\"+str(int(i))+r\"}$\"  for i in ticks])\n",
    "            \n",
    "            axL=plt.axes([bbax3.x0, bbax0.y0+bbax0.height*0.1,  bbax3.width, bbax0.height*.5], frameon=False, xticks=[],yticks=[])\n",
    "            axL.set_xlim((0,1))\n",
    "            axL.set_ylim((0,1))\n",
    "            ellipse = patches.Ellipse(xy=(0.1,0.8), height=0.2, width=0.15, angle=0 ,      edgecolor='red', fc='None', lw=1)\n",
    "            axL.text(0.35, 0.75, \"1-sigma CR\")\n",
    "            axL.add_patch(ellipse)\n",
    "            axC=plt.axes([bbax3.x0, bbax0.y0+bbax0.height*0.1,bbax3.width,bbax0.height*0.1], frameon=True)\n",
    "            norm = mpl.colors.Normalize(vmin=0,vmax=1)\n",
    "            sm = plt.cm.ScalarMappable(cmap=plt.cm.viridis, norm=norm)\n",
    "            sm.set_array([])\n",
    "            cbar = plt.colorbar(sm,cax=axC,orientation='horizontal')\n",
    "            cbar.set_ticks([0,1])\n",
    "            cbar.set_ticklabels(['low','high'])\n",
    "            axC.set_title('Probability density')\n",
    "            axC.xaxis.get_major_ticks()[0].label1.set_horizontalalignment('left')\n",
    "            axC.xaxis.get_major_ticks()[1].label1.set_horizontalalignment('right')\n",
    "            #f.suptitle('hepato ' + n,fontsize=30)\n",
    "            plt.tight_layout()    \n",
    "        f.savefig(fimag,dpi=600)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "considerable-negotiation",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "million-gossip",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "induced-gabriel",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "thorough-mount",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "graphic-trial",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "instructional-egyptian",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "subsequent-suicide",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "biological-given",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "superior-stability",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "present-buddy",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "opposite-intro",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "compressed-delicious",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "c14",
   "language": "python",
   "name": "c14"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}