diff --git a/tutorials/HFF_100KB/.README.md.swp b/tutorials/HFF_100KB/.README.md.swp new file mode 100644 index 0000000..ae0e3bf Binary files /dev/null and b/tutorials/HFF_100KB/.README.md.swp differ diff --git a/tutorials/HFF_100KB/NuclearDeformation/.gitignore b/tutorials/HFF_100KB/NuclearDeformation/.gitignore index 05e5ce4..7cace18 100644 --- a/tutorials/HFF_100KB/NuclearDeformation/.gitignore +++ b/tutorials/HFF_100KB/NuclearDeformation/.gitignore @@ -1 +1,2 @@ simulation.py +.ipynb_checkpoints diff --git a/tutorials/HFF_100KB/README.md b/tutorials/HFF_100KB/README.md index d4ad010..4d7e3f9 100644 --- a/tutorials/HFF_100KB/README.md +++ b/tutorials/HFF_100KB/README.md @@ -1,6 +1,6 @@ ## HFF model with nuclear landmarks at 100KB resolution -This folder provides tutorials for setting up and performing simulations of the HFF nucleusi at the 100KB resolution with the presence of various nuclear landmarks, including nuclear lamina, speckles, and nucleoli. +This folder provides tutorials for setting up and performing simulations of the HFF nucleus at the 100KB resolution with the presence of various nuclear landmarks, including nuclear lamina, speckles, and nucleoli. - SphericalNucleus 1. This tutorial explains how to setup simulations of a spherical nucleus. diff --git a/tutorials/HFF_100KB/SphericalNucleus/.gitignore b/tutorials/HFF_100KB/SphericalNucleus/.gitignore index 05e5ce4..7cace18 100644 --- a/tutorials/HFF_100KB/SphericalNucleus/.gitignore +++ b/tutorials/HFF_100KB/SphericalNucleus/.gitignore @@ -1 +1,2 @@ simulation.py +.ipynb_checkpoints diff --git a/tutorials/HFF_100KB/SphericalNucleus/.ipynb_checkpoints/simulation-checkpoint.ipynb b/tutorials/HFF_100KB/SphericalNucleus/.ipynb_checkpoints/simulation-checkpoint.ipynb deleted file mode 100644 index a659b4d..0000000 --- a/tutorials/HFF_100KB/SphericalNucleus/.ipynb_checkpoints/simulation-checkpoint.ipynb +++ /dev/null @@ -1,190 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "abbcbe8f-b163-46a4-9c8f-84d194be4309", - "metadata": {}, - "source": [ - "## Import the necessary packages" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cce78b79-5e3c-4712-8cd8-2fd3875c2cb0", - "metadata": {}, - "outputs": [], - "source": [ - "import parmed as pmd\n", - "import json\n", - "import sys\n", - "from sys import platform\n", - "import mdtraj as md\n", - "import mdtraj.reporters\n", - "import random\n", - "from openNucleome import OpenNucleome" - ] - }, - { - "cell_type": "markdown", - "id": "0987436b-c228-4ba7-baa9-28e21b3a63ed", - "metadata": {}, - "source": [ - "## Important parameters\n", - "We defined the important parameters in the next block. First, we set the transition probability between dP particles and P particles as 0.2, and the transition frequency as 4000. When creating the system, we set type 6, 7 as the dP and P particles, respectively, so we kept using 6, 7 here. In this example, we ran a simulation with total length of 3,000,000 steps, and output one configuration and the energy every 2000 steps." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "157264b7-9a55-48d5-b369-b389aa292dd8", - "metadata": {}, - "outputs": [], - "source": [ - "prob_P_dP = 0.2 # Transition probability from P to dP\n", - "prob_dP_P = 0.2 # Transition probability from dP to P\n", - "transition_freq = 4000\n", - "sampling_freq = 2000\n", - "dP_type = 6\n", - "P_type = 7\n", - "total_steps = 3000000" - ] - }, - { - "cell_type": "markdown", - "id": "8e0e9ef8-97a2-488e-ba82-5fdfe446b9b3", - "metadata": {}, - "source": [ - "## Initialize the system\n", - "We first set up an example \"model\" of class \"OpenNucleome\" with the conserved temperature, damping coefficient, timestep, and the mass scale. In this folder, we also included the initial configuration \"human.pdb\" used for the simulation and created a system according to the initial configuration.\n", - "\n", - "Because during the simulations, we would transit the dP speckle and P speckle particles, here, we logged the start index and end index of speckle particle, and computed the number of speckle particles. \n", - "\n", - "In this example, we freezed all the lamina beads, and would not consider the dynamics of the membrane, and that is why we set \"False\" (off) for membrane dynamics. Consequently, there was also no need to set the bond between lamina beads, so we set \"None\" for the variable \"membrane_bond\"." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8ba4e26b-cacb-48ec-93c9-e6fa6409a9d1", - "metadata": {}, - "outputs": [], - "source": [ - "model = OpenNucleome(1.0, 0.1, 0.005, 1.0) # 1.0: temperature (LJ reduced unit); \n", - " # 0.1: damping coefficient (LJ reduced unit);\n", - " # 0.005: timestep (LJ reduced unit);\n", - " # 1.0: mass_scale\n", - " \n", - "PDB_file = \"human.pdb\" #The initial configuration\n", - "start_spec_index = model.N_chr_nuc+1\n", - "end_spec_index = model.N_chr_nuc_spec+1\n", - "N_spec = start_spec_index-end_spec_index #Currently, the number of speckle is 1600.\n", - "\n", - "# Generate new elements and construct topology as well\n", - "# membrane_dynamics: True for including lamina dynamics, False for excluding lamina dynamics;\n", - "# membrane_bond: A file contains the lamina bond when membrane_dynamics is on.\n", - "model.create_system(PDB_file, membrane_dynamics = False, membrane_bond = None) " - ] - }, - { - "cell_type": "markdown", - "id": "07c89f5c-7839-45ae-8b98-186d5600b42e", - "metadata": {}, - "source": [ - "## Add the force field\n", - "In this example, we loaded default settings for the interactions between chromosomes and chromosomes, and chromosomes and nuclear landmarks. All the types of potential can be found in \"Section: Energy function of the whole nucleus model\" in Supporting Information. According to the order of added potential, the index of speckle-speckle potential is 6, and this index is needed when transiting the speckle particle between type 6 and type 7." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "44d8e300-07d3-4e37-8991-9820bffb5a58", - "metadata": {}, - "outputs": [], - "source": [ - "# Add the default force field\n", - "model.load_default_settings()\n", - "\n", - "index_spec_spec_potential = 6" - ] - }, - { - "cell_type": "markdown", - "id": "6c942879-33b9-4ddd-92bf-0ebae07561cb", - "metadata": {}, - "source": [ - "## Perform the simulation\n", - "We first created the simulation with a specific Platform, and here, we used \"CUDA\" but users can also use \"CPU\", \"Reference\", and \"OpenCL\" according to their hardware. Before the simulation, we minimized the energy to make the system much more reasonable and stable. After randomly setting velocity, we started our simulation with a total length of 3,000,000 steps and output the configuration and energy every 2000 steps, and change the speckle types every 4000 steps as we mentioned in the previous blocks" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5f55005e-3048-44ea-9a53-cbc4a55e2940", - "metadata": {}, - "outputs": [], - "source": [ - "#model.save_system(\"model_before_simulation_0.xml\")\n", - "\n", - "simulation = model.create_simulation(platform_type = \"CUDA\") # Users can also use CPU, Reference, OpenCL.\n", - "simulation.context.setPositions(model.chr_positions)\n", - "\n", - "simulation.minimizeEnergy()\n", - "\n", - "simulation.reporters.append(mdtraj.reporters.DCDReporter('HFF_3e6_every2000.dcd', sampling_freq))\n", - "\n", - "def setVelocity(context):\n", - " sigma = u.sqrt(1.0*u.kilojoule_per_mole / model.chr_system.getParticleMass(1))\n", - " velocs = u.Quantity(1.0 * np.random.normal(size=(model.chr_system.getNumParticles(), 3)),\n", - " u.meter) * (sigma / u.meter)\n", - " context.setVelocities(velocs)\n", - "setVelocity(simulation.context)\n", - "\n", - "simulation.reporters.append(mmapp.statedatareporter.StateDataReporter(sys.stdout, sampling_freq, step=True,\n", - " potentialEnergy=True, kineticEnergy=True, temperature=True, progress=True,\n", - " remainingTime=True, separator='\\t', totalSteps = total_steps))\n", - "\n", - "for i in range(total_steps//transition_freq):\n", - " simulation.step(transition_freq)\n", - " # Change the type of speckles every 4000 steps, non-equilibrium scheme.\n", - "\n", - " # Do the chemical modification, and change the spec-spec potential on the fly.\n", - " for j in np.random.randint(start_spec_index, end_spec_index, N_spec): \n", - "\n", - " if model.compart_type[j] == dP_type-1:\n", - " model.compart_type[j] = P_type-1 if random.random() < prob_dP_P else dP_type-1\n", - " else:\n", - " model.compart_type[j] = dP_type-1 if random.random() < prob_P_dP else P_type-1\n", - "\n", - " # Update the context after changing the type of speckles.\n", - " for m in range(model.chr_system.getNumParticles()):\n", - " model.chr_system.getForce(index_spec_spec_potential).setParticleParameters(m, [model.compart_type[m]])\n", - " model.chr_system.getForce(index_spec_spec_potential).updateParametersInContext(simulation.context)\n", - "\n", - "# Keep the final result of spec types in case constructing the configuration for the continuous simulation.\n", - "np.savetxt('compt_final_frame.txt', (np.array(model.compart_type)+1).reshape((-1,1)), fmt='%d')" - ] - } - ], - "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.8" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/tutorials/HFF_100KB/SphericalNucleus/README.md b/tutorials/HFF_100KB/SphericalNucleus/README.md new file mode 100644 index 0000000..0b08207 --- /dev/null +++ b/tutorials/HFF_100KB/SphericalNucleus/README.md @@ -0,0 +1,5 @@ +## Spherical Nucleus + +This folder provides tutorials for setting up and performing simulations of the HFF nucleus at the 100KB resolution with the presence of various nuclear landmarks, including nuclear lamina, speckles, and nucleoli. + + diff --git a/tutorials/analysis_code/README.md b/tutorials/analysis_code/README.md index a31dd33..6722ebb 100644 --- a/tutorials/analysis_code/README.md +++ b/tutorials/analysis_code/README.md @@ -1,13 +1,13 @@ ## Analysis code -compute_contact_prob.ipynb: call and computes the contact probabilities (DamID, TSASeq, and chr-chr) +- compute_contact_prob.ipynb: call and computes the contact probabilities (DamID, TSASeq, and chr-chr) -contact_calculation: execution file compiled from Fortran code "../../openNucleome/utility/chromosome_contact_calculation.f90" +- contact_calculation: execution file compiled from Fortran code "../../openNucleome/utility/chromosome_contact_calculation.f90" -mol_info: store all the molecular information files necessary for computing contact probabilities +- mol_info: store all the molecular information files necessary for computing contact probabilities -frame_10.dcd: trajectory file example which contains 10 frames +- frame_10.dcd: trajectory file example which contains 10 frames -contact_prob: store all the chromosome-chromosome contact probabilities, the number of contacts, and the number of analyzed frames +- contact_prob: store all the chromosome-chromosome contact probabilities, the number of contacts, and the number of analyzed frames -expt_constraints_HFF_100KB.txt: experimental chromosome-chromosome contact probabilities, the first 2489 values are ideal contacts, the next 6 are compt-compt contacts, and the remaining 231 are interchromosomal contacts +- expt_constraints_HFF_100KB.txt: experimental chromosome-chromosome contact probabilities, the first 2489 values are ideal contacts, the next 6 are compt-compt contacts, and the remaining 231 are interchromosomal contacts diff --git a/tutorials/analysis_code/contact_prob/README.md b/tutorials/analysis_code/contact_prob/README.md index f08bcb3..3fcb822 100644 --- a/tutorials/analysis_code/contact_prob/README.md +++ b/tutorials/analysis_code/contact_prob/README.md @@ -1,7 +1,7 @@ ## contact_prob -contact_prob.txt: averaged contact probabilities over different frames +- contact_prob.txt: averaged contact probabilities over different frames -nframes.txt: the number of analyzed frames +- nframes.txt: the number of analyzed frames -counter.txt: the total number of contacts +- counter.txt: the total number of contacts diff --git a/tutorials/analysis_code/mol_info/README.md b/tutorials/analysis_code/mol_info/README.md index 42427bb..716c6e7 100644 --- a/tutorials/analysis_code/mol_info/README.md +++ b/tutorials/analysis_code/mol_info/README.md @@ -1,17 +1,17 @@ ## Molecular information -DAM-ID_HFF_100KB.txt: DamID signal for HFF model at resolution of 100KB +- DAM-ID_HFF_100KB.txt: DamID signal for HFF model at resolution of 100KB -TSA-Seq_SON_HFF_100KB.txt: TSA-Seq signal for HFF model at resolution of 100KB +- TSA-Seq_SON_HFF_100KB.txt: TSA-Seq signal for HFF model at resolution of 100KB -gLengthFile.txt: The difference between each two neighboring values in the file represents the length of the chromosome +- gLengthFile.txt: The difference between each two neighboring values in the file represents the length of the chromosome -maternalIdxFile.txt: Index of each maternal chromosome +- maternalIdxFile.txt: Index of each maternal chromosome -paternalIdxFile.txt: Index of each paternal chromosome +- paternalIdxFile.txt: Index of each paternal chromosome -compartment_genome_100KB_diploid_HFF.txt: The type of each chromatin bead (1: compartment A, 2: compartment B, 3: compartment C, 4: beads around type 3) +- compartment_genome_100KB_diploid_HFF.txt: The type of each chromatin bead (1: compartment A, 2: compartment B, 3: compartment C, 4: beads around type 3) -mol_index_cell_type-HFF.txt: Index of the chromosome that each bead belongs to +- mol_index_cell_type-HFF.txt: Index of the chromosome that each bead belongs to -tad_index_genome_100KB_diploid_cell_type-HFF.txt: The first column is the index of each bead in the diploid representation and the second column is the index in the haploid representation +- tad_index_genome_100KB_diploid_cell_type-HFF.txt: The first column is the index of each bead in the diploid representation and the second column is the index in the haploid representation