non-hepatocyte_and_combined.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "occupied-client",
   "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": "specific-conspiracy",
   "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": "genetic-constitutional",
   "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": "potential-noise",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "consecutive-welding",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('C14data_liver_samples_20211020.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "humanitarian-employment",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Data files can be found in Zenodo \n",
    "\n",
    "directory_uns = 'liver_uns/'\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 == \"unsorted\" and pathology != [\"Y\", \"C\"]')\n",
    "data['sigma14C'] = data['e14C']\n",
    "data_uns = data\n",
    "\n",
    "\n",
    "directory_non = 'liver_non/'\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['Dbirth'] -= 0.5\n",
    "data['age'] = data['Dcoll'] - data['Dbirth']\n",
    "#data = data.query('type == \"hepatocyte\" and pathology != \"Y\"')\n",
    "data = data.query('type == \"non-hepatocyte\" and pathology != [\"Y\",\"C\"]')\n",
    "data['sigma14C'] = data['e14C']\n",
    "data_non = data\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "reflected-distinction",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_non"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "imposed-chaos",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.seterr(all='ignore')\n",
    "results= dict()\n",
    "NUM=0\n",
    "results= dict(data=data_uns,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']:\n",
    "        continue\n",
    "    path = directory_uns+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",
    "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_unsorted.xlsx')\n",
    "calc_c14_all(results)\n",
    "printt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "personalized-quality",
   "metadata": {},
   "outputs": [],
   "source": [
    " rename={'sigma':r'$\\sigma$','d':r'$\\delta$'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "irish-technology",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "for n,v in results['results'].items():  \n",
    "    fimag = 'figs/unsorted'+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",
    "                        logparas=results['models'][n].logparas)    \n",
    "        f = plt.gcf()\n",
    "        f.show();\n",
    "        #f.suptitle('hepato ' + n,fontsize=30)\n",
    "        f.savefig(fimag,dpi=600)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bacterial-relative",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.seterr(all='ignore')\n",
    "NUM=0\n",
    "results_non= dict(data=data_non,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']:\n",
    "        continue\n",
    "    path = directory_non+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_non['results'].update(  {m.__class__.__name__+'spline_' : {'raw':res} } )\n",
    "    results_non['models'].update(   {m.__class__.__name__+'spline_' : m } )        \n",
    "\n",
    "get_arviz(results_non,burnin=1000)\n",
    "ranking_non = get_getranking(results_non)\n",
    "ranking_non\n",
    "printt_non = get_pointestimates(results_non,error_in_real=False)[['Parameter','Value','Confidence Interval']]\n",
    "printt_non.to_excel('liver_non.xlsx')\n",
    "calc_c14_all(results_non)\n",
    "printt_non"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "directed-permit",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "for n,v in results_non['results'].items():  \n",
    "    fimag = 'figs/non'+n+'.png'\n",
    "    if not os.path.isfile(fimag) or True:\n",
    "        c14.corner.corner_R(v['azdata'].posterior.a.values,\n",
    "                   results_non['models'][n].parameter_names,\n",
    "                        point_estimate=results_non['models'][n].transform_physical_to_fit(v['median'],mode='bayes'),\n",
    "                        rename=rename,\n",
    "                        logparas=results_non['models'][n].logparas)    \n",
    "        f = plt.gcf()\n",
    "        f.show();\n",
    "        #f.suptitle('hepato ' + n,fontsize=30)\n",
    "        f.savefig(fimag,dpi=600)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "convenient-portrait",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "impressed-symbol",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "\n",
    "subject = 85\n",
    "excel_writer = pd.ExcelWriter('time_rates_flows_uns.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",
    "    df = pd.DataFrame(ipid)\n",
    "    df.index.name = 'age'\n",
    "    df.index =df.index/10\n",
    "    \n",
    "    ipid.update({n:np.array([v]) for n,v in parameter_phy.items()})\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",
    "excel_writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "located-health",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "particular-trading",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "\n",
    "subject = 68\n",
    "excel_writer = pd.ExcelWriter('time_rates_flows_non.xlsx')\n",
    "for m_i,(model_name,model) in enumerate(results_non['models'].items()):\n",
    "    print(model_name)\n",
    "    age = results_non['data'].loc[subject,'age']\n",
    "    T = np.arange(0,90,0.1)\n",
    "    parameter_phy =  {i:results_non['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_non['data'].index) if x == subject][0]  \n",
    "    ipid = listofdict_to_dictofarray_f(ipsd,a_i)\n",
    "    df = pd.DataFrame(ipid)\n",
    "    df.index.name = 'age'\n",
    "    df.index =df.index/10\n",
    "    \n",
    "    ipid.update({n:np.array([v]) for n,v in parameter_phy.items()})\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",
    "excel_writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "equal-missouri",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "devoted-gibraltar",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "civil-stroke",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "recovered-resistance",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "equivalent-occasions",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "under-championship",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "czech-opportunity",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "other-transmission",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "modern-black",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "behavioral-encoding",
   "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
}