Single Chain: Chromatin Dynamics Simulations on Chromosome 10 of GM12878 Cell Line

This tutorial should take between 5 to 15 minutes of reading and performing simulations.

The first step is to import the OpenMiChroM module

from OpenMiChroM.ChromDynamics import MiChroM
from OpenMiChroM.CndbTools import cndbTools

MiChroM class sets the initial parameters of the simulation:

time_step=0.01(Simulation time step to perfom the integration) temperature=1.0 (Set the temperature of your simulation)

sim = MiChroM(temperature=1.0, time_step=0.01)

There are three hardware platform options to run the simulations:

platform="cuda"
platform="opencl"
platform="cpu"

sim.setup(platform="opencl") # Double-check CUDA installation in your system

Sets the directory name where to save the simulation outputs

sim.saveFolder('output_chr10')

The next step is to load the chromatin type annotations sequence for chromosome 10 and generate an initial 3D structure to start the simulation.

The first column should contain the locus index. The second column should have the locus type annotation. A template file of the chromatin sequence of types can be found here.

Here we set the chromosome polymer initial 3D configuration using the create_springSpiral function.

The loci positions are stored in the variable Chrom10 as a NumPy array \([N:3]\), where \(N\) is the number of beads.

Chrom10 = sim.create_springSpiral(ChromSeq='inputs/chr10_beads.txt', isRing=False)

print(Chrom10[:5]) #Print the firsts five beads positions.

Load the initial structure into the sim object.

The center=True moves your system to the origin.

sim.loadStructure(Chrom10, center=True)

The chromosome 3D initial structure can be saved in .ndb file format. The file is stored in the path given in saveFolder.

sim.saveStructure(mode = 'ndb')

The next step is to add the force field in the simulation object sim.

In this tutorial, the forces are added in two sets:

MiChroM Homopolymer (Bonded) Potentials

sim.addFENEBonds(kfb=30.0)
sim.addAngles(ka=2.0)
sim.addRepulsiveSoftCore(Ecut=4.0)

MiChroM non-Bonded Potentials

sim.addTypetoType(mu=3.22, rc = 1.78)
sim.addIdealChromosome(mu=3.22, rc = 1.78, dinit=3, dend=500)

The last potential adds a spherical restraint to collapse the initial structure.

sim.addFlatBottomHarmonic( kr=5*10**-3, n_rad=15.0)

Run a short simulation to generate a collapsed structure.

block = 3*10**2
n_blocks = 2*10**3

Two variables control the chromatin dynamics simulations steps:

block: The number of steps performed in each cycle (n_blocks) n_blocks: The number of blocks of the simulation.

The initial collapse simulation will run for \(3\times10^2 \times 2\times10^3 = 6\times10^5\) steps.

for _ in range(n_blocks):
    sim.runSimBlock(block, increment=False)

Details about the output of each simulation block:

• bl=0 The index number of the simulated block. The parameter increment=False, is used to ignore the steps counting.

• pos[1]=[X,Y,Z] The spatial position for the locus \(1\).

• dr=1.26 The average of the loci displacements in each block (in units os sigma).

• t=0 The current time step of the simulated. The parameter increment=False ignores the steps counting and keeps \(t=0.0\).

• kin=1.5 The kinetic energy of the system.

• pot=19.90 The total potential energy of the system.

• RG=7.654 The Radius of gyration at the end of the simulated block.

• SPS=12312 The Steps Per Second of each block.

The radius of gyration is a good parameter to check if the collapse was well performed. If the chromosome polymer is not collapsed, it is necessary to rerun the initial collapse steps.

print(sim.chromRG())
sim.saveStructure(mode='ndb')
#sim.saveStructure((mode='gro') #There are options to save the chromosome 3D structure using the standard file format used for macromolecules.
#sim.saveStructure(mode='pdb')

The next step is to remove the spherical restraint force to run the production simulation. The OpenMM commands are: sim.system.getForces() Returns a list of all forces in sim system. sim.system.removeForce(5) Removes the force indexed as 5 (the last one in this tutorial).

sim.system.getForces()

sim.system.removeForce(5)

To run the production simulation, it is necessary to initialize the .cndb file to save the chromatin dynamics trajectory.

sim.initStorage(filename="traj_chr10")

Sets the parameters of the production simulations:

\(block = 5\times10^2\) \(n\_blocks = 2\times10^3\)

block = 5*10**2
n_blocks = 2*10**3

for _ in range(n_blocks):
    sim.runSimBlock(block, increment=True)
    sim.saveStructure()

Once the simulation is performed, it is necessary to close the .cndb file to avoid losing the trajectory data.

sim.storage[0].close()

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).

Trajectory Analyses using cndbTools

cndbTools is a class that allows analyses in the chromatin dynamics trajectories using the binary format .cndb (compact ndb).

cndbTools = cndbTools()

Load the cndb file in the variable chr10_traj.

chr10_traj = cndbTools.load('output_chr10/traj_chr10_0.cndb')

print(chr10_traj) # Print the information of the cndb trajectory.

Extract the loci XYZ position over the simulated 2000 frames and save in the variable chr10_xyz.

chr10_xyz = cndbTools.xyz(frames=[1,2000,1], beadSelection='all', XYZ=[0,1,2])

The variable chr10_xyz allows the cndbTools to perform several analyses. In this example, the Radius of Gyration can be obtained as a function of the simulated frames.

import matplotlib.pyplot as plt
import matplotlib as mpl


chr10_RG = cndbTools.compute_RG(chr10_xyz)
plt.plot(chr10_RG)
plt.ylabel(r'Radius of Gyration ($\sigma$)',fontsize=15)
plt.xlabel(r'Simulation Frames',fontsize=15)

cndbTools allows the selection of beads to compute the analyses. An example is the Radial Distribution Probability (RDP) for each chromatin type (A and B).

chr10_A = cndbTools.xyz(frames=[1,2000,1], beadSelection=chr10_traj.dictChromSeq['A1'], XYZ=[0,1,2])
chr10_B = cndbTools.xyz(frames=[1,2000,1], beadSelection=chr10_traj.dictChromSeq['B1'], XYZ=[0,1,2])

print("Computing RDP...")
r_A, RDP_chr10_A = cndbTools.compute_RDP(chr10_A, radius=15.0, bins=200)
r_B, RDP_chr10_B = cndbTools.compute_RDP(chr10_B, radius=15.0, bins=200)

plt.plot(r_A, RDP_chr10_A, color='red', label='A')
plt.plot(r_B, RDP_chr10_B, color='blue', label='B')
plt.xlabel(r'r ($\sigma$)', fontsize=15,fontweight='normal', color='k')
plt.ylabel(r'$\rho(r)/N_{type}$', fontsize=15,fontweight='normal', color='k')
plt.legend()
plt.gca().set_xlim([1/200,15.0])

cndbTools also generates the in silico Hi-C map (contact probability matrix).

In this tutorial, the trajectory contains 2000 snapshots of chromosome 10 of the GM12878 cell line. It is expected the in silico Hi-C not be fully converged due to the inadequate sampling. To produce a converged map, it is recommended to simulate around 20 replicas with 10000 frames on each, which generates an ensemble of 200.000 chromosome structures.

print("Generating the contact probability matrix...")
chr10_sim_HiC = cndbTools.traj2HiC(chr10_xyz)

plt.matshow(chr10_sim_HiC, norm=mpl.colors.LogNorm(vmin=0.001, vmax=chr10_sim_HiC.max()),cmap="Reds")
plt.colorbar()

To visualize the chromosome’s 3D structures in the standard visualization software 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.

The ndb plugin for visualizing the chromatin dynamics trajectories in VMD/Chimera/Pymol is under development.