3_bbmodel_concordance_v2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compared to v1, v2 also calculates the 95% two-sided CI for rho, by basic bootstrap method with 100 iterations (can be very slow, as the MLE optimisation is repeated 100 times) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/giulio/Desktop/Generational_multiplex_analysis/\n"
     ]
    }
   ],
   "source": [
    "import itertools as itt\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import scipy.optimize as spo\n",
    "import scipy.special as spsp\n",
    "import scipy.stats as sps\n",
    "import seaborn as sns\n",
    "import time\n",
    "\n",
    "from matplotlib.ticker import FixedLocator\n",
    "\n",
    "mpl.rcParams['axes.titlesize'] = 'xx-large'\n",
    "mpl.rcParams['axes.labelsize'] = 'xx-large'\n",
    "mpl.rcParams['xtick.labelsize'] = 'x-large'\n",
    "mpl.rcParams['ytick.labelsize'] = 'x-large'\n",
    "mpl.rcParams['xtick.direction'] = 'out'\n",
    "mpl.rcParams['ytick.direction'] = 'out'\n",
    "mpl.rcParams['legend.frameon'] = True\n",
    "mpl.rcParams['legend.framealpha'] = 0.5\n",
    "mpl.rcParams['legend.fontsize'] = 'large'\n",
    "\n",
    "path_proj = os.getcwd()+ '/'\n",
    "path_sc = path_proj + 'csv/Single_cell/'\n",
    "path_sort = path_proj + 'csv/Sort/'\n",
    "path_plot = path_proj + 'figures/'\n",
    "print(path_proj)\n",
    "\n",
    "#LOAD THE PROCESSED DATA\n",
    "df = pd.read_csv(path_sc+'Pooled_data.csv', sep=';', decimal=',')\n",
    "\n",
    "#GLOBAL VARIABLES FOR ITERATION\n",
    "conds = ['P', 'P+ILs']\n",
    "or_cells = ['SLAM-HSC', 'ST-HSC', 'MPP']\n",
    "\n",
    "#SET A RANDOM SEED TO REPLICATE THE BOOTSTRAP CONFIDENCE INTEVALS\n",
    "#BASED ON RANDOM SAMPLING WITH REPLACEMENT\n",
    "np.random.seed(123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#IMPORTANT: SAMPLING PROBABILITY MUST BE SPECIFIED FROM THE USER! \n",
    "r_sample = 0.71"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#NOTE: A FULL FAMILY VECTOR IS THE VECTOR OF CELL COUNTS PER GENERATION (ONE GENERATION PER VECTOR ENTRY);\n",
    "#A SAMPLED FAMILY VECTOR IS A FULL FAMILY VECTOR WHERE CELLS WERE NOT NECESSARILY COUNTED/OBSERVED\n",
    "\n",
    "#BELOW ARE THE FUNCTIONS AND CLASSES TO DEFINE FULL AND SAMPLED FAMILY VECTORS\n",
    "\n",
    "#GIVEN A SAMPLED FAMILY VECTOR, POPULATE THE NEXT GENERATION WITH ALL POSSIBLE COUNTS THAT ARE PHYSICALLY POSSIBLE\n",
    "#(HERE IT IS ASSUMED THAT EACH DIVIDING CELLS PRODUCES TWO CELLS IN THE NEXT GENERATION)\n",
    "def add_generation(family):\n",
    "    i = len(family)\n",
    "    past = [family[j] * 2**(i-j) for j in range(i)]\n",
    "    max_val = 2**i - sum(past)\n",
    "    new_fam_list = [family[:]+[k] for k in range(max_val+1)]\n",
    "    return new_fam_list\n",
    "\n",
    "#CREATE ALL SAMPLED FAMILY VECTORS UP TO GENERATION DD\n",
    "def fam_generator(DD):\n",
    "    possible_fam = [[0],[1]]   \n",
    "    for DD in range(DD):\n",
    "        temp = [add_generation(elem) for elem in possible_fam]\n",
    "        possible_fam = list(itt.chain.from_iterable(temp))\n",
    "    return possible_fam\n",
    "    \n",
    "#GIVEN A FULL FAMILY VECTOR, RETURNS ALL SAMPLED FAMILY VECTORS THAT MAY POSSIBLY DESCENT FROM IT  \n",
    "def sub_trees_gen(ft):\n",
    "    res_lst = [[float(k)] for k in range(int(ft[0])+1)]\n",
    "    for max_n in ft[1:]:\n",
    "        new_lst = [vec[:]+[float(k)] for k in range(int(max_n)+1)\n",
    "                                     for vec in res_lst]\n",
    "        res_lst = new_lst\n",
    "    return [Sub_tree(np.asarray(vec), sum(np.subtract(ft,vec)))\n",
    "                for vec in res_lst if sum(vec)!=0.]\n",
    "\n",
    "#CLASS FOR SAMPLED FAMILY VECTORS\n",
    "class Sub_tree(object):\n",
    "    def __init__(self, st, lost_c):\n",
    "        self.st = st\n",
    "        self.lost_c = lost_c\n",
    "        self.range = np.arange(len(st))[st>0].max() - np.arange(len(st))[st>0].min()\n",
    "    def show_st(self):\n",
    "        print( 'ST:', self.st, 'Lc:', self.lost_c)\n",
    "\n",
    "#CLASS FOR FULL FAMILY VECTORS\n",
    "class Full_tree(object):\n",
    "    def __init__(self, ft):\n",
    "        self.ft = ft\n",
    "        self.ebc = cells_arrived_in_gen(ft)\n",
    "        self.st_lst = sub_trees_gen(ft)\n",
    "        self.range = np.arange(len(ft))[ft>0].max() - np.arange(len(ft))[ft>0].min()\n",
    "    def show(self):\n",
    "        print( 'FT:', self.ft)\n",
    "        print( 'Ebc:', self.ebc)\n",
    "        for st in self.st_lst:\n",
    "            st.show_st()\n",
    "\n",
    "            \n",
    "#BETA-BINOMIAL (BB) MODEL FOR CELL DIVISION,\n",
    "#WHERE A FAMILY DEVELOPS FROM A 1 CELL IN GENERATION 0 ITERATIVELY AS FOLLOWS:\n",
    "#THE NUMBER OF CELLS THAT DIVIDE FROM GENERATION i TO i+1 FOLLOWS A BETA-BINOMIAL DISTRIBUTION\n",
    "#PARAMETRISED SO THAT ONE CELL DIVIDES WITH PROBABILITY p, AND THE CELLS IN THE SAME GENERATION\n",
    "#DECIDE TO DIVIDE WITH (INTRACLASS) CORRELATION rho.\n",
    "#IN ACCORDANCE WITH EXPERIMENTAL DATA, THE CELLS FROM FULL FAMILY VECTOR GROWN UNDER THE BB MODEL\n",
    "#ARE THEN SAMPLED (INDEPENDENTLY AND WITH EQUAL DISTRIBUTION) WITH PROBABILITY r.\n",
    "\n",
    "#GIVEN A FULL FAMILY VECTOR m, RETURN THE NUMBER OF CELLS THAT WERE EVER BORN IN EACH GENERATION\n",
    "def cells_arrived_in_gen(m):\n",
    "    nn = np.zeros(len(m))\n",
    "    nn[0] = 1.\n",
    "    for k in range(1, len(m)):\n",
    "        nn[k] = 2.*(nn[k-1] - m[k-1])\n",
    "    return nn\n",
    "\n",
    "#PROBABILITY THAT A BETA-BINOMIAL DISTRIBUTION, WITH PARAMETERS n, a, b, EQUALS k\n",
    "def bb_dist(k, n, a, b):\n",
    "    if a == 0.:\n",
    "        if k==0.:\n",
    "            return 1.\n",
    "        else:\n",
    "            return 0.\n",
    "    elif b == 0.:\n",
    "        if k == n:\n",
    "            return 1.\n",
    "        else:\n",
    "            return 0.\n",
    "    else:\n",
    "        return spsp.binom(n,k) * spsp.beta(k+a,n-k+b) / spsp.beta(a,b)\n",
    "    \n",
    "#PROBABILITY THAT A BINOMIAL DISTRIBUTION, WITH PARAMETERS n, R, EQUALS k    \n",
    "def bino_dist(k, n, r):\n",
    "    if r == 0.:\n",
    "        if k==0.:\n",
    "            return 1.\n",
    "        else:\n",
    "            return 0.\n",
    "    elif r == 1.:\n",
    "        if k == n:\n",
    "            return 1.\n",
    "        else:\n",
    "            return 0.\n",
    "    else:\n",
    "        return spsp.binom(n, k) * r**(k) * (1.-r)**(n-k)\n",
    "    \n",
    "#PROBEBILITY OF OBSERVING THE FULL FAMILY VECTOR m UNDER THE BB MODEL WITH PARAMETERS p, rho\n",
    "def pmf_m(m, p, rho):\n",
    "    nn = cells_arrived_in_gen(m)\n",
    "    if rho==0.0:\n",
    "        res = np.prod([bino_dist(nn[i]-m[i], nn[i], p[i]) for i in range(len(m)-1)])\n",
    "        return res\n",
    "    elif rho==1.:\n",
    "        m_v = [(m[dd] / (2.**dd)) == 1. for dd in range(len(m))]\n",
    "        if sum(m_v) == 1.:\n",
    "            n_dd = next(dd for dd in range(len(m)) if (m[dd]/(2.**dd))==1.)\n",
    "            res = np.prod([p[i] for i in range(n_dd)])*(1.-p[n_dd])\n",
    "            return res\n",
    "        else:\n",
    "            return 0.\n",
    "    else:    \n",
    "        alpha = np.multiply(p, (1.-rho)/rho)\n",
    "        beta = np.multiply(np.subtract(1., p), (1.-rho)/rho)    \n",
    "        res = np.prod([bb_dist(m[i], nn[i], beta[i], alpha[i])\n",
    "                       for i in range(len(alpha)-1)])\n",
    "        return res\n",
    "\n",
    "#PROBABILITY OF OBSERVING THE SAMPLED FAMILY VECTOR v,\n",
    "#SAMPLING EACH CELL FROM THE FULL FAMILY VECTOR m WITH SUCCESS RATE r\n",
    "def pmf_v_given_m(v, m, r):\n",
    "    return np.prod([bino_dist(v[i], m[i], r)\n",
    "                    for i in range(len(v))]) / (1.-((1.-r)**sum(m)))\n",
    "\n",
    "#PROBABILITY OF OBSERVING FRAY/RANGE EQUAL TO k UNDER THE BB MODEL rho, p,\n",
    "#GENERATING FULL FAMILY VECTORS IN Ft_lst, THEN SAMPLED WITH SUCCESS RATE r \n",
    "def fray_dist(k, rho, p, r, Ft_lst):\n",
    "    out = 0.\n",
    "    for Ft in Ft_lst:\n",
    "        out += sum([pmf_v_given_m(St.st, Ft.ft, r) for St in Ft.st_lst\n",
    "                    if St.range==k]) * pmf_m(Ft.ft, p, rho)\n",
    "    return out\n",
    "\n",
    "#NEGATIVE LOG-LIKELIHOOD OF FRAY/RANGE DATA\n",
    "def log_like_data(rho, p, r, Ft_lst, fray_freq):\n",
    "    return -sum([fray_freq[k]*np.log(fray_dist(k, rho, p, r, Ft_lst))\n",
    "                 for k in range(len(fray_freq))])\n",
    "\n",
    "#BELOW ARE THE FUNCTIONS TO CALCLUATE THE DISTRIBUTION AND CONFIDENCE INTERVALS\n",
    "#OF DATA IN A GIVEN SUPPORT (SIMILARLY AS IN NOTEBOOK 2)\n",
    "\n",
    "#TRANSFORM THE PANDAS DATAFRAME OF A FAMILY INTO A FAMILY VECTOR\n",
    "def df_fam_to_vec(dff, max_g):\n",
    "    res = np.zeros(max_g+1)\n",
    "    gens, counts = np.unique(dff.Generation, return_counts=True)\n",
    "    res[gens] = counts\n",
    "    return res\n",
    "\n",
    "def epmf_on_support(data, support, bool_norm=True):\n",
    "    classes, counts = np.unique(data, return_counts=True)\n",
    "    counts_sum = counts.sum()\n",
    "    if counts_sum == 0:\n",
    "        return np.zeros(len(support))\n",
    "    elif bool_norm:\n",
    "        return np.hstack([counts[classes==cl] if any(classes==cl) else np.zeros(1) for cl in support]\n",
    "                        )/counts_sum\n",
    "    else:\n",
    "        return np.hstack([counts[classes==cl] if any(classes==cl) else np.zeros(1) for cl in support])\n",
    "    \n",
    "def stat_bootCI(func, args, data, n_iter=250000, bool_trim01=True):\n",
    "    stat = func(data, *args)\n",
    "    boot_stat = np.array([func(data[boot_index], *args)\n",
    "                          for boot_index in np.random.choice(len(data), size=(n_iter,len(data)))])\n",
    "    boot_per = np.percentile(boot_stat, q=[97.5, 2.5], axis=0)\n",
    "    boot_ci = 2.*stat - boot_per\n",
    "    if bool_trim01:\n",
    "        boot_ci[0,boot_ci[0,:]<0.] = np.zeros((boot_ci[0,:]<0.).sum())\n",
    "        boot_ci[1,boot_ci[1,:]>1.] = np.ones((boot_ci[1,:]>1.).sum())\n",
    "    return stat, boot_ci"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " ST-HSC P 77\n",
      "Full trees generation time: 0.017768122999996194\n",
      "\n",
      "Optimization terminated successfully;\n",
      "The returned value satisfies the termination criteria\n",
      "(using xtol =  1e-05 )\n",
      "MLE rho estimation time: 0.4749418570000046\n",
      "LR test pval= 9.876525720395996e-12\n",
      "Best-fit rho: {} in ({},{})\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'NoneType' object has no attribute 'format'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-11-8c8fd264fac4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     80\u001b[0m         \u001b[0mboot_per_rho\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpercentile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mboot_rho\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m97.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2.5\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     81\u001b[0m         \u001b[0mboot_ci_rho\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m2.\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mest_rho\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mboot_per_rho\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 82\u001b[0;31m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Best-fit rho: {} in ({},{})'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mround\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mest_rho\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mround\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mboot_ci_rho\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mround\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mboot_ci_rho\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     83\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     84\u001b[0m         \u001b[0;31m#PLOTTING CODE\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'NoneType' object has no attribute 'format'"
     ]
    }
   ],
   "source": [
    "###FIG 2F\n",
    "df = df[df.Culture_time=='48h']#Use 48h data only\n",
    "boot_iter = 250000\n",
    "for n_oc,oc in enumerate(or_cells):\n",
    "    for n_cnd,cnd in enumerate(conds):\n",
    "        #if n_oc != 1 or n_cnd != 0:\n",
    "        #    continue\n",
    "            \n",
    "        df_temp = df[(df.Original_cell==oc)&(df.Culture_condition==cnd)]\n",
    "        \n",
    "        #SET MAX GENERATION AND ALL SAMPLED FAMILY VECTORS FROM THE DATA \n",
    "        max_g = df_temp.Generation.max()\n",
    "        uni_fams = np.unique(df_temp.Family)\n",
    "        print('\\n', oc, cnd, len(uni_fams))\n",
    "        fam_vec = np.array([df_fam_to_vec(df_temp[df_temp.Family==fam], max_g)\n",
    "                            for fam in uni_fams])\n",
    "        \n",
    "        #ESTIMATE p FOR THE BB MODEL, THE PROBABILITY FOR A CELL TO DIVIDE TO NEXT GENERATION \n",
    "        cc_pop = fam_vec.sum(axis=0)*(2.** -np.arange(max_g+1))\n",
    "        p_div = np.array([cc_pop[i+1:].sum()/cc_pop[i:].sum() for i in np.arange(len(cc_pop)-1)] + [0.])\n",
    "        \n",
    "        #CALCULATE THE OBSERVED FRAY/RANGE, AND ITS DISTRIBUTION, FROM THE DATA\n",
    "        fray_data = np.array([df_temp[df_temp.Family==fam].Generation.max()\n",
    "                              - df_temp[df_temp.Family==fam].Generation.min()\n",
    "                              for fam in uni_fams])\n",
    "        fray_support = np.arange(max_g)\n",
    "        fray_freq = epmf_on_support(fray_data, fray_support, bool_norm=False)\n",
    "        \n",
    "        #COMMENT/UNCOMMENT TO RECALCULATE THE OPTIMISATION STEP\n",
    "        fray_epmf, boot_ci = stat_bootCI(func=epmf_on_support, args=(fray_support,), data=fray_data)\n",
    "\n",
    "        #Generate all full amily vectors up to max observed depth d\n",
    "        #start = time.process_time()\n",
    "        samp_vec = np.array([fam for fam in fam_generator(max_g) if (fam * 2.**-np.arange(len(fam))).sum()<=1.\n",
    "                                                    and (fam * 2.**-np.arange(len(fam))).sum()>=0.])\n",
    "        full_vec = np.array([fam for fam in samp_vec if (fam * 2.**-np.arange(len(fam))).sum()==1.])\n",
    "        Ft_lst = [Full_tree(fam) for fam in full_vec]\n",
    "        #end = time.process_time()\n",
    "        #print('Full trees generation time:', end-start)\n",
    "        \n",
    "        #MLE rho estimation\n",
    "        #start = time.process_time()\n",
    "        solution = spo.minimize_scalar(log_like_data, bounds=(0., 1.),\n",
    "                            args=(p_div, r_sample, Ft_lst, fray_freq),\n",
    "                            method='bounded', options={'disp': True})\n",
    "        est_rho = solution.x\n",
    "        #end = time.process_time()\n",
    "        #print('MLE rho estimation time:', end-start)\n",
    "        \n",
    "        #BOOTSTRAP RHO\n",
    "        boot_n_rho = 100\n",
    "        boot_rho = []\n",
    "        for boot_data in np.random.choice(fray_data, size=(boot_n_rho,len(fray_data))):\n",
    "            boot_freq = epmf_on_support(boot_data, fray_support, bool_norm=False)\n",
    "            boot_solution = spo.minimize_scalar(log_like_data, bounds=(0., 1.),\n",
    "                                    args=(p_div, r_sample, Ft_lst, boot_freq),\n",
    "                                    method='bounded', options={'disp': False})\n",
    "            boot_rho.append(boot_solution.x)\n",
    "    \n",
    "        #Fray distribution under the BB model for other values of rho\n",
    "        rho_lst = np.array([est_rho, 0.9, 0.5, 0.0])\n",
    "        fray_pmf = {str(rho):[fray_dist(gen, rho, p_div, r_sample, Ft_lst)\n",
    "                              for gen in fray_support] for rho in rho_lst}\n",
    "        \n",
    "        #LIKELIHOOD-RATIO TEST (EXTRA)\n",
    "        neg_log_like_H0 = log_like_data(0., p_div, r_sample, Ft_lst, fray_freq)\n",
    "        neg_log_like_H1 = log_like_data(est_rho, p_div, r_sample, Ft_lst, fray_freq)\n",
    "        like_ratio_stat = 2. * (neg_log_like_H0 - neg_log_like_H1)\n",
    "        pval = sps.chi2.sf(like_ratio_stat, df=1)\n",
    "        print('LR test pval=', pval)\n",
    "        \n",
    "        #SAVE VARIABLES THAT TAKE A LONG TIME TO CALCULATE\n",
    "        with open('./pickled_data/Fray_ci_'+oc+cnd, 'wb') as fp:\n",
    "            pickle.dump((fray_epmf,boot_ci,est_rho,boot_rho,fray_pmf,pval), fp)\n",
    "        #COMMENT/UNCOMMENT ENDS HERE\n",
    "        \n",
    "        with open('./pickled_data/Fray_ci_'+oc+cnd, 'rb') as fp:\n",
    "            fray_epmf,boot_ci,est_rho,boot_rho,fray_pmf,pval = pickle.load(fp)\n",
    "\n",
    "        #PRINT RHO CI\n",
    "        boot_per_rho = np.percentile(boot_rho, q=[97.5, 2.5], axis=0)\n",
    "        boot_ci_rho = 2.*est_rho - boot_per_rho\n",
    "        print('Best-fit rho: {} in 95% CI ({},{})'.format(round(est_rho,3),round(boot_ci_rho[0],3),round(boot_ci_rho[1],3)))            \n",
    "        \n",
    "        #PLOTTING CODE\n",
    "        rho_lst = np.array([est_rho, 0.9, 0.5, 0.0])###THIS IS ADDED FOR WHEN THE PART ABOVE IS COMMENTED\n",
    "        width = 0.5/len(rho_lst)\n",
    "        t_cols = sns.color_palette('Blues', n_colors=6)[::-1]\n",
    "        color_dct = {str(est_rho):t_cols[0], '0.9':t_cols[1],\n",
    "                     '0.5':t_cols[2], '0.0':t_cols[3]}\n",
    "        shift_dct = {str(rho_lst[k]):-0.30 + (k+1)*0.125 for k in range(len(rho_lst))}    \n",
    "        sns.set_style(\"ticks\")\n",
    "        fig, ax = plt.subplots()\n",
    "        ax.bar(x=fray_support-0.30, height=fray_epmf, width=width, align='edge',\n",
    "           color='Red', label='Data')\n",
    "        for rho in rho_lst:\n",
    "            if rho == est_rho:\n",
    "                ax.bar(x=fray_support+shift_dct[str(rho)], height=fray_pmf[str(rho)],\n",
    "                       width=width, align='edge', color=color_dct[str(rho)], hatch='///',\n",
    "                       label=r'$\\rho = $'+str(round(rho, 3))+' (Best fit)')\n",
    "            else:\n",
    "                ax.bar(x=fray_support+shift_dct[str(rho)], height=fray_pmf[str(rho)],\n",
    "                       width=width, align='edge', color=color_dct[str(rho)],\n",
    "                       label=r'$\\rho = $'+str(round(rho, 3)))     \n",
    "        xx_fray = np.add(np.arange(len(fray_epmf)), (width/2.)-0.30)\n",
    "        for k in range(len(fray_epmf)):\n",
    "            ax.plot([xx_fray[k]]*2, [boot_ci[0,k], boot_ci[1,k]],\n",
    "                     color='Black', ls='-', lw=1., mfc='Black', mec='Black',\n",
    "                     marker='_', mew=1, ms=5.)\n",
    "        sns.despine(fig=fig)\n",
    "        ax.set_xlabel('Range', fontsize='xx-large')\n",
    "        ax.set_ylabel('% of clones', fontsize='xx-large')\n",
    "        ax.set_xticks(np.arange(max_g))#len(fray_epmf)))\n",
    "        ax.set_xticklabels(np.arange(max_g), fontsize='x-large')# np.arange(len(fray_epmf))\n",
    "        ax.set_yticks(np.arange(0., 1.1, 0.25))\n",
    "        ax.set_yticklabels(('0', '25', '50', '75', '100'), fontsize='x-large')\n",
    "        if oc == 'SLAM-HSC' and cnd == 'P':\n",
    "            ax.set_xlim(-0.5, 2.5)\n",
    "        elif oc == 'SLAM-HSC' and cnd == 'P+ILs':\n",
    "            ax.set_xlim(-0.5, 3.5)\n",
    "        elif oc == 'ST-HSC' and cnd == 'P':\n",
    "            ax.set_xlim(-0.5, 3.5)\n",
    "        elif oc == 'ST-HSC' and cnd == 'P+ILs':\n",
    "            ax.set_xlim(-0.5, 4.5)\n",
    "        elif oc == 'MPP' and cnd == 'P':\n",
    "            ax.set_xlim(-0.5, 2.5)\n",
    "        elif oc == 'MPP' and cnd == 'P+ILs':\n",
    "            ax.set_xlim(-0.5, 4.5)\n",
    "        #ax.set_xlim(-0.5, fray_data.max()+0.5)\n",
    "        ax.set_ylim(bottom=-0.0, top=1.0)\n",
    "        ax.legend(fontsize='x-large', loc='upper left', ncol=1, bbox_to_anchor= (1.05, 1.))\n",
    "        ax.set_title('Generation range 48h '+oc+' '+cnd, fontdict = {'fontsize':'xx-large'})\n",
    "        fig.savefig(path_plot+'/Fray fit 48h '+oc+' '+cnd+'.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best-fit rho: 0.799 in 95% CI (0.668,0.947)\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}