{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Single Chromosome Simulation" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "This tutorial should take between 20 to 30 minutes of reading and performing simulations." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "#### Chromatin Dynamics Simulations on Chromosome 10 of GM12878 Cell Line" ] }, { "attachments": {}, "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" ] }, { "attachments": {}, "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
\n", "- `temperature=1.0`: set the temperature of your simulation
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sim = MiChroM(temperature=1.0, timeStep=0.01)" ] }, { "attachments": {}, "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\")\n", "sim.setup(platform=\"cuda\")" ] }, { "attachments": {}, "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')" ] }, { "attachments": {}, "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/inputs/chr10_beads.txt).
\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)" ] }, { "attachments": {}, "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])" ] }, { "attachments": {}, "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)" ] }, { "attachments": {}, "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)" ] }, { "attachments": {}, "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)" ] }, { "attachments": {}, "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, nRad=15.0)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Run a short simulation to generate a collapsed structure." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Creates the simulation context and initializes the system\n", "sim.createSimulation()" ] }, { "attachments": {}, "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": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "block = 3*10**2\n", "n_blocks = 2*10**3" ] }, { "attachments": {}, "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)
\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": {}, "outputs": [], "source": [ "for step in range(n_blocks):\n", " sim.run(block)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Details about the output of each simulation block:\n", "\n", "- `Step=0`: index number of the simulated block.
\n", "- `RG=7.654`: radius of gyration at the end of the simulated block.
\n", "- `Etotal=19.90`: total energy of the system (reduced units).
\n", "- `Epot=19.90`: total potential energy of the system (reduced units).
\n", "- `Ekin=1.5`: kinetic energy of the system (reduced units).
" ] }, { "attachments": {}, "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": [ "sim.saveStructure(mode='ndb')" ] }, { "attachments": {}, "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')" ] }, { "attachments": {}, "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()" ] }, { "attachments": {}, "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\")" ] }, { "attachments": {}, "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.createReporters(statistics=True, traj=True, trajFormat=\"cndb\", energyComponents=True, interval=5*10**2)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Set the parameters of the production simulation:\n", "\n", "$block = 5\\times10^2$
\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 }, "outputs": [], "source": [ "for step in range(n_blocks):\n", " sim.run(block)" ] }, { "attachments": {}, "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)." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "#### Trajectory analysis using cndbTools" ] }, { "attachments": {}, "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()" ] }, { "attachments": {}, "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/OpenMiChroM_0.cndb')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(chr10_traj) # Print the information of the cndb trajectory." ] }, { "attachments": {}, "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=range(0,2000,1), XYZ=[0,1,2])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "max([int(key) for key in chr10_traj.cndb.keys() if key != 'types'])" ] }, { "attachments": {}, "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)" ] }, { "attachments": {}, "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=range(0,2000,1), beadSelection=chr10_traj.dictChromSeq[b'A1'], XYZ=[0,1,2])\n", "chr10_B1 = cndbTools.xyz(frames=range(0,2000,1), beadSelection=chr10_traj.dictChromSeq[b'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])" ] }, { "attachments": {}, "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()" ] }, { "attachments": {}, "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": "OpenMichrom Latest", "language": "python", "name": "openmm2" }, "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.12.0" }, "vscode": { "interpreter": { "hash": "cf0b1156a9c1fd50aa79dcca714d7a34c3b0f1748625c1528654ed7c64ee5692" } } }, "nbformat": 4, "nbformat_minor": 4 }