01_example_for_running_random_boundary_simulations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/hbrandao/anaconda3/lib/python3.7/site-packages/Cython/Compiler/Main.py:369: FutureWarning: Cython directive 'language_level' not set, using 2 for now (Py2). This will change in a later release! File: /home/hbrandao/libs/MCMs-as-random-barriers-to-loop-extrusion-paper/entropic_smcTranslocator.pyx\n",
      "  tree = Parsing.p_module(s, pxd, full_module_name)\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.stats import expon\n",
    "\n",
    "import pyximport;\n",
    "pyximport.install()\n",
    "\n",
    "from entropic_smcTranslocator import smcTranslocatorDirectional \n",
    "\n",
    "from brandaolib.contactProbability_generator import ChainLayout\n",
    "from itertools import product\n",
    "\n",
    "\n",
    "import pandas as pd\n",
    "import os\n",
    "import re\n",
    "from scipy.sparse import coo_matrix\n",
    "\n",
    "from multiprocessing import Pool\n",
    "import pickle\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define loop strength calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getLoopStrength(Z,Zcount,loopBases,pad=16,shift=50):\n",
    "    ZZ = Z/Zcount\n",
    "    ZZ[np.isnan(ZZ)] = 0\n",
    "    np.fill_diagonal(ZZ,np.mean(np.diag(ZZ))/2)\n",
    "    data = ZZ + ZZ.T\n",
    "    \n",
    "    # skip first and last TADs\n",
    "    loopStrengthsList = [] \n",
    "    for (stBin,endBin) in zip(loopBases[2:],loopBases[1:-1]):\n",
    "        MM = data[stBin - pad:stBin + pad, endBin - pad:endBin + pad]\n",
    "        MC = data[stBin - pad+shift:stBin + pad+shift, endBin - pad+shift:endBin + pad+shift]\n",
    "        \n",
    "        M = MM/MC\n",
    "        \n",
    "        M[np.isinf(M)] = 0\n",
    "        # divide box into 3 equal spaces\n",
    "        L = len(M)\n",
    "        box1 = M[0:L//3,0:L//3] \n",
    "        box2 = M[L//2-L//6:L//2+L//6,L//2-L//6:L//2+L//6] \n",
    "        box3 = M[L-L//3:L,L-L//3:L] \n",
    "        \n",
    "        loopStrengthsList.append(np.nansum(box2)/(np.nansum(box3)+np.nansum(box1))*2)\n",
    "    return loopStrengthsList\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define range of parameter sweep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# fixed boundaries\n",
    "nreps = 1000 #10000 #20000\n",
    "nsamples = 70\n",
    "fractionBoundaries = [0]#[2,4,8,16]#[0,0.24,0.5,1,1.6,2]\n",
    "boundaryPauseProb = [0]#[0.5,0.25,0.125,0.0625,0.03125]#[1,0.75,0.5,0.25]\n",
    "processivities = [60,75,95,100]#[100,150,200,250,300]#[100,200,300]\n",
    "separations = [75, 100, 150,200]#[100,150,200,250,300]#[200]\n",
    "\n",
    "TAD_size = 300\n",
    "TAD_boundaries = np.arange(0,10000,TAD_size)\n",
    "\n",
    "name_format_string = './RandomBoundarySims_refined/FixedTADs_proc{}_sep{}_tadsize{}_fractionBoundaries{}_boundaryPauseProb{}_nsamples{}_nreps{}.pkl'\n",
    "lof_nums = list(product(fractionBoundaries,boundaryPauseProb,processivities,separations))\n",
    "lof = [(x[0],x[1],name_format_string.format(x[2],x[3],TAD_size,x[0],x[1],nreps,nsamples),x[2],x[3]) for x in lof_nums]\n",
    "\n",
    "\n",
    "len(lof)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # random boundaries\n",
    "# nreps = 20000\n",
    "# nsamples = 70\n",
    "# fractionBoundaries = [4,8,16]#[0,0.24,0.5,1,1.6,2]\n",
    "# boundaryPauseProb = [0.5,0.25,0.125,0.0625]#[1,0.75,0.5,0.25]\n",
    "\n",
    "# lof_nums = list(product(fractionBoundaries,boundaryPauseProb))\n",
    "\n",
    "# lof = [(x[0],x[1],'./RandomBoundarySims/fractionBoundaries{}_boundaryPauseProb{}_nsamples{}_nreps{}.pkl'\n",
    "#                                 .format(x[0],x[1],nreps,nsamples)) for x in lof_nums]\n",
    "\n",
    "\n",
    "# TAD_size = 200\n",
    "# TAD_boundaries =##np.random.choice(10000,int(10000/TAD_size),replace=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the simulation for all parameters (parallel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc60_sep200_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc75_sep200_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc95_sep200_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc100_sep200_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc60_sep150_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc75_sep150_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc100_sep150_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc95_sep150_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc75_sep100_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc95_sep100_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc100_sep100_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc60_sep100_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc100_sep75_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc95_sep75_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc60_sep75_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n",
      "Doing ./RandomBoundarySims_refined/FixedTADs_proc75_sep75_tadsize300_fractionBoundaries0_boundaryPauseProb0_nsamples1000_nreps70.pkl\n"
     ]
    }
   ],
   "source": [
    "def doOne(idx):\n",
    "    try:\n",
    "        filename = lof[idx][2]\n",
    "\n",
    "        #unchanging parameters\n",
    "        v_extruder =  50 # kb/min  \n",
    "        BELT_ON=0\n",
    "        BELT_OFF=1\n",
    "        switchRate = 0 \n",
    "        SWITCH_PROB= switchRate # switching rate\n",
    "        PUSH=0\n",
    "        PAIRED=0\n",
    "        SLIDE=1\n",
    "        SLIDE_PAUSEPROB=0.99 # what is this value?\n",
    "        loop_prefactor=1.5 # what is this value?\n",
    "        FULL_LOOP_ENTROPY=1 # what is this value?\n",
    "        FRACTION_ONESIDED=0\n",
    "\n",
    "        # Extruder dynamics frequency of sampling parameters #\n",
    "        numExtruderSteps = 1000 # steps taken for each simulation sample\n",
    "        numInitializationSteps = 10000 # how long we take to equilibrate the simulation\n",
    "\n",
    "\n",
    "        # Polymer and extruder dynamics parameters #\n",
    "        L = 10000 \n",
    "        processivity = lof[idx][3]#100\n",
    "        separations = lof[idx][4]#2*processivity\n",
    "#         TAD_size = 200\n",
    "#         TAD_boundaries = np.arange(0,L,TAD_size)\n",
    "        PAUSEPROB=0.0 # motor pause probability\n",
    "\n",
    "\n",
    "        smcNum = L//separations # number of SMCs loaded  \n",
    "        SWITCH =  np.ones(L,dtype=np.double)*SWITCH_PROB\n",
    "        LIFETIME = processivity\n",
    "        birthArray = np.ones(L)/L\n",
    "        deathArray = np.zeros(L, dtype=np.double) + 1. / LIFETIME\n",
    "        deathArray[0:1] = 1\n",
    "        deathArray[L-2:L-1] = 1\n",
    "        stallDeathArray = deathArray \n",
    "        stallDeathArray[0:1] = 1\n",
    "        stallDeathArray[L-2:L-1] = 1    \n",
    "        pauseArray = np.zeros(L, dtype=np.double) + PAUSEPROB \n",
    "        slidePauseArray = np.zeros(L, dtype=np.double) + SLIDE_PAUSEPROB\n",
    "        oneSidedArray = np.zeros(smcNum, dtype=np.int64)\n",
    "        belt_on_array = np.zeros(smcNum, dtype=np.double) + BELT_ON\n",
    "        belt_off_array = np.zeros(smcNum, dtype=np.double) + BELT_OFF\n",
    "        spf=slidePauseArray*(1.-(1.-SLIDE_PAUSEPROB)*np.exp(-1.*loop_prefactor))\n",
    "        spb=slidePauseArray*(1.-(1.-SLIDE_PAUSEPROB)*np.exp(loop_prefactor))\n",
    "        ################### TAD BOUNDARY CONDITION###################\n",
    "        stallLeftArray = np.zeros(L, dtype = np.double)\n",
    "        stallRightArray = np.zeros(L,dtype = np.double)\n",
    "        stallprob = 0.4\n",
    "        for b in TAD_boundaries:\n",
    "            stallLeftArray[b] = stallprob\n",
    "            stallRightArray[b] = stallprob\n",
    "        ##################################################################\n",
    "        ################### Random barrier CONDITION###################\n",
    "        numRandomBarriers = int(np.round(len(TAD_boundaries)*lof[idx][0]))\n",
    "        randomBarriers = np.random.choice(L,numRandomBarriers,replace=False)\n",
    "        boundaryPauseProb = lof[idx][1]\n",
    "        for b in randomBarriers:\n",
    "            stallLeftArray[b] = boundaryPauseProb\n",
    "            stallRightArray[b] = boundaryPauseProb\n",
    "\n",
    "        ##################################################################    \n",
    "\n",
    "\n",
    "        transloc = smcTranslocatorDirectional(birthArray, deathArray, stallLeftArray, stallRightArray, pauseArray,\n",
    "                                             stallDeathArray, smcNum, oneSidedArray, FRACTION_ONESIDED, slide=SLIDE,\n",
    "                                             slidepauseForward=spf, slidepauseBackward=spb, switch=SWITCH, pushing=PUSH,\n",
    "                                            belt_on=belt_on_array, belt_off=belt_off_array,SLIDE_PAUSEPROB=SLIDE_PAUSEPROB) \n",
    "\n",
    "        transloc.steps(numInitializationSteps)\n",
    "\n",
    "\n",
    "        ### dump empty data as test\n",
    "        X = []\n",
    "        Y = []\n",
    "        P = []\n",
    "        C = []\n",
    "        loopSizes_array = []\n",
    "        L_tot = L\n",
    "        Z = coo_matrix((P,(X,Y)),shape=(L_tot,L_tot))\n",
    "        Zcount = coo_matrix((C,(X,Y)),shape=(L_tot,L_tot))\n",
    "        ChIP = np.zeros(L)\n",
    "\n",
    "        param_sweep_dict = {'ChIP':ChIP,'Z':Z,'Zcount':Zcount,\\\n",
    "                    'RandomBarriers':randomBarriers,\\\n",
    "                   'nSMCs':smcNum,\\\n",
    "                   'numReps':nreps,'nsamples_per_map':nsamples,'lifetime':LIFETIME,\\\n",
    "                        'TAD_boundaries':TAD_boundaries,\\\n",
    "                       'randomBoundaryPauseProb':boundaryPauseProb,'TADBoundaryPauseProb':stallprob,'loopSizesArray':loopSizes_array}\n",
    "\n",
    "\n",
    "        pickle.dump(param_sweep_dict,open(filename,'wb'))\n",
    "\n",
    "        print(\"Doing {}\".format(filename))\n",
    "\n",
    "\n",
    "\n",
    "        for ic in range(nreps):\n",
    "            transloc.steps(numExtruderSteps)\n",
    "            # generate chain using loop and gap statistics\n",
    "            smcs_lr = transloc.getSMCs()\n",
    "            \n",
    "            # get ChIP-seq\n",
    "            for x in range(len(smcs_lr[0])):                \n",
    "                ChIP[smcs_lr[0][x]] += 1\n",
    "                ChIP[smcs_lr[1][x]] += 1\n",
    "            \n",
    "            \n",
    "            smcs = [(smcs_lr[0][x], smcs_lr[1][x]) for x in range(len(smcs_lr[0]) )  ]\n",
    "            cl = None\n",
    "            try:\n",
    "                if len(smcs)==0:\n",
    "                    smcs = [(0,1)]\n",
    "                cl = ChainLayout(smcs,L_tot)\n",
    "            except:\n",
    "                print(\"Pseudoknot formed count:{}\".format(ic))\n",
    "                print(cl)\n",
    "                assert(1==0)\n",
    "                cl = None\n",
    "            if cl is None:\n",
    "                continue\n",
    "                \n",
    "            # measure loop sizes\n",
    "            loopSizes_array = loopSizes_array + list(smcs_lr[1]-smcs_lr[0])\n",
    "            \n",
    "            # get contact maps, subsampling from distribution \n",
    "            vals = sorted(np.random.choice(L_tot,nsamples, replace=False))\n",
    "            for ix in range(len(vals)):\n",
    "                for iy in range(ix,len(vals)):\n",
    "                    x = vals[ix]\n",
    "                    y = vals[iy]\n",
    "                    deff = cl.get_dist(x,y)\n",
    "                    if deff == 0:\n",
    "                        pc = 1\n",
    "                    else:\n",
    "                        pc = 1/np.sqrt(deff)**3\n",
    "                    if not np.isnan(pc):\n",
    "                        X.append(x)\n",
    "                        Y.append(y) # x\n",
    "                        P.append(pc) # probability\n",
    "                        C.append(1) # counts\n",
    "\n",
    "        Z = coo_matrix((P,(X,Y)),shape=(L_tot,L_tot))\n",
    "        Zcount = coo_matrix((C,(X,Y)),shape=(L_tot,L_tot))\n",
    "\n",
    "\n",
    "        param_sweep_dict = {'ChIP':ChIP,'Z':Z,'Zcount':Zcount,\\\n",
    "                    'RandomBarriers':randomBarriers,\\\n",
    "                   'nSMCs':smcNum,\\\n",
    "                   'numReps':nreps,'nsamples_per_map':nsamples,'lifetime':LIFETIME,\\\n",
    "                        'TAD_boundaries':TAD_boundaries,\\\n",
    "                       'randomBoundaryPauseProb':boundaryPauseProb,'TADBoundaryPauseProb':stallprob,'loopSizesArray':loopSizes_array}\n",
    "\n",
    "        pickle.dump(param_sweep_dict,open(filename,'wb'))\n",
    "        \n",
    "        return 1\n",
    "    except:\n",
    "        return 0\n",
    "#     return 1\n",
    "\n",
    "pool = Pool(len(lof))\n",
    "pool.map(doOne, np.arange(len(lof)))\n",
    "# pool = Pool(30)\n",
    "# pool.map(doOne, np.arange(len(lof)))"
   ]
  },
  {
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}