docs/source/Tutorials/Tutorial_Single_Chromosome.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single Chromosome Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial should take between 20 to 30 minutes of reading and performing simulations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Chromatin Dynamics Simulations on Chromosome 10 of GM12878 Cell Line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to import the **OpenMiChroM** module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from OpenMiChroM.ChromDynamics import MiChroM\n",
    "from OpenMiChroM.CndbTools import cndbTools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MiChroM` class sets the initial parameters of the simulation:\n",
    "\n",
    "- `time_step=0.01`: set the simulation time step to perfom the integration<br>\n",
    "- `temperature=1.0`: set the temperature of your simulation<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = MiChroM(temperature=1.0, time_step=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are four hardware platform options to run the simulations: \n",
    "```python\n",
    "platform=\"cuda\"\n",
    "platform=\"opencl\"\n",
    "platform=\"hip\"\n",
    "platform=\"cpu\"\n",
    "```\n",
    "\n",
    "Choose accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.setup(platform=\"opencl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the directory name in which the output of the simulation is saved:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.saveFolder('output_chr10')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to load the chromatin compartment sequence for chromosome 10 and generate an initial 3D structure to start the simulation. We can use the [createSpringSpiral](https://open-michrom.readthedocs.io/en/latest/OpenMiChroM.html#OpenMiChroM.ChromDynamics.MiChroM.createSpringSpiral) function to set the initial configuration of the polymer based in the sequence file.\n",
    "\n",
    "The first column of the sequence file should contain the locus index. The second should have the locus type annotation. A template file of the chromatin sequence of types can be found [here](https://github.com/junioreif/OpenMiChroM/blob/main/Tutorials/Chromosomes_simulations/inputs/chr10_beads.txt).<br>\n",
    "\n",
    "The loci positions are stored in the variable **chr10** as a NumPy array $[N:3]$, where $N$ is the number of beads. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chr10 = sim.createSpringSpiral(ChromSeq='inputs/chr10_beads.txt', isRing=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the position of the first five beads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(chr10[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial structure should then be loaded into the `sim` object.\n",
    "\n",
    "The option `center=True` moves your system to the origin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.loadStructure(chr10, center=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial 3D chromosome structure can be saved in [.ndb file format](https://ndb.rice.edu/ndb-format). The file is stored in the path given in `saveFolder`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.saveStructure(mode='ndb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to add the force field in the simulation object `sim`.\n",
    "\n",
    "In this tutorial, the forces can be divided into two sets:\n",
    "\n",
    "**MiChroM Homopolymer (Bonded) Potentials**  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.addFENEBonds(kfb=30.0)\n",
    "sim.addAngles(ka=2.0)\n",
    "sim.addRepulsiveSoftCore(Ecut=4.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**MiChroM Non-Bonded Potentials**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.addTypetoType(mu=3.22, rc=1.78)\n",
    "sim.addIdealChromosome(mu=3.22, rc=1.78, dinit=3, dend=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last potential adds a spherical constrain to collapse the initial structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.addFlatBottomHarmonic(kr=5*10**-3, n_rad=15.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run a short simulation to generate a collapsed structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "block = 3*10**2\n",
    "n_blocks = 2*10**3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two variables control the chromatin dynamics simulation steps:\n",
    "\n",
    "`block`: The number of time-steps performed in each cycle (or block)</br>\n",
    "`n_blocks`: The number of cycles (or blocks) simulated. \n",
    "\n",
    "The initial collapse simulation will run for  $3\\times10^2 \\times  2\\times10^3 = 6\\times10^5$ time-steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for _ in range(n_blocks):\n",
    "    sim.runSimBlock(block, increment=False) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Details about the output of each simulation block:\n",
    "\n",
    "- `bl=0`: index number of the simulated block. The parameter `increment=False` is used to ignore the steps counting.<br>\n",
    "- `pos[1]=[X,Y,Z]`: spatial position for the locus 1. <br>\n",
    "- `dr=1.26`: average of the loci displacements in each block (in units of sigma). <br>\n",
    "- `t=0`: current simulation time. <br>\n",
    "- `kin=1.5`: kinetic energy of the system (reduced units). <br>\n",
    "- `pot=19.90`: total potential energy of the system (reduced units). <br>\n",
    "- `RG=7.654`: radius of gyration at the end of the simulated block. <br>\n",
    "- `SPS=12312`: steps per second of each block."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The radius of gyration is a good parameter to check the performance of the collapse.\n",
    "If the chromosome polymer is not collapsed, it is necessary to rerun the initial collapse steps. We can also save the structure for inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sim.chromRG())\n",
    "sim.saveStructure(mode='ndb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The structure can also be saved using stardard file formats used for macromolecules, as the `pdb` and `gro` formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.saveStructure(mode='gro')\n",
    "sim.saveStructure(mode='pdb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to remove the spherical constrain force to run the production simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.removeFlatBottomHarmonic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If necessary, one could remove any of the forces applied in the system. To see the forces in the system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.forceDict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sim.removeForce(forceName=\"TypetoType\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run the production simulation, it is necessary to initialize the .cndb file to save the chromatin dynamics trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.initStorage(filename=\"traj_chr10\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the parameters of the production simulation:\n",
    "\n",
    "$block = 5\\times10^2$ </br> \n",
    "$n\\_blocks = 2\\times10^3$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "block = 5*10**2\n",
    "n_blocks = 2*10**3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "for _ in range(n_blocks):\n",
    "    sim.runSimBlock(block, increment=True) \n",
    "    sim.saveStructure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the simulation is completed, it is necessary to close the .cndb file to avoid losing the trajectory data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.storage[0].close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation should generate the `traj_chr10_0.cndb` trajectory file in the output_chr10 folder. This file contains 2000 frames (one snapshot per block)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Trajectory analysis using cndbTools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`cndbTools` is a class that allows analyses in the chromatin dynamics trajectories using the binary format [.cndb](https://ndb.rice.edu/ndb-format) (compact ndb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cndbTools = cndbTools()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the cndb file in the variable `chr10_traj`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chr10_traj = cndbTools.load('output_chr10/traj_chr10_0.cndb')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(chr10_traj) # Print the information of the cndb trajectory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the loci XYZ position over the simulated 2000 frames and save in the variable `chr10_xyz`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chr10_xyz = cndbTools.xyz(frames=[1,2001,1], beadSelection='all', XYZ=[0,1,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable `chr10_xyz` allows the cndbTools to perform several analyses.\n",
    "In this example, the radius of gyration can be obtained as a function of the simulated frames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "\n",
    "chr10_RG = cndbTools.compute_RG(chr10_xyz)\n",
    "plt.plot(chr10_RG)\n",
    "plt.ylabel(r'Radius of Gyration ($\\sigma$)',fontsize=11)\n",
    "plt.xlabel(r'Simulation Frames',fontsize=11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`cndbTools` allows the selection of beads to compute the analyses. An example is the Radial Distribution Probability (RDP) for each chromatin subcompartments A1 and B1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chr10_A1 = cndbTools.xyz(frames=[1,2001,1], beadSelection=chr10_traj.dictChromSeq['A1'], XYZ=[0,1,2])\n",
    "chr10_B1 = cndbTools.xyz(frames=[1,2001,1], beadSelection=chr10_traj.dictChromSeq['B1'], XYZ=[0,1,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Computing RDP...\")\n",
    "r_A1, RDP_chr10_A1 = cndbTools.compute_RDP(chr10_A1, radius=15.0, bins=200)\n",
    "r_B1, RDP_chr10_B1 = cndbTools.compute_RDP(chr10_B1, radius=15.0, bins=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(r_A1, RDP_chr10_A1, color='red', label='A')\n",
    "plt.plot(r_B1, RDP_chr10_B1, color='blue', label='B')\n",
    "plt.xlabel(r'r ($\\sigma$)', fontsize=11,fontweight='normal', color='k')\n",
    "plt.ylabel(r'$\\rho(r)/N_{type}$', fontsize=11,fontweight='normal', color='k')\n",
    "plt.legend()\n",
    "plt.gca().set_xlim([1/200,15.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use `cndbTools` to generate the *in silico* Hi-C map (contact probability matrix).\n",
    "\n",
    "In this tutorial, the trajectory contains 2,000 snapshots of chromosome 10 of the GM12878 cell line. For this set of structures, we expect the *in silico* Hi-C to not be fully converged due to inadequate sampling. \n",
    "To produce a converged map, it is recommended to simulate around 20 replicas with 10,000 frames on each, which generates an ensemble of 200,000 chromosome structures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Generating the contact probability matrix...\")\n",
    "chr10_sim_HiC = cndbTools.traj2HiC(chr10_xyz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.matshow(chr10_sim_HiC, norm=mpl.colors.LogNorm(vmin=0.001, vmax=chr10_sim_HiC.max()),cmap=\"Reds\")  \n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize the chromosome's 3D structures in the standard visualization softwares for macromolecules, there are available scripts for converting the `ndb`/`cndb` file format to `.pdb` and `.gro`. For details, please check the [Nucleome Data Bank](https://ndb.rice.edu/ndb-format).\n",
    "\n",
    "The `ndb` plugin for visualizing the chromatin dynamics trajectories in VMD/Chimera/Pymol is under development."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.11 ('base')",
   "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.11"
  },
  "vscode": {
   "interpreter": {
    "hash": "2b315b238527074b48ede1db88856294695fedbd9a1befcadf019011d6f4a552"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}