add new features in utils

openNucleome/README.md

@@ -12,4 +12,4 @@
 
 - lamina.py: All the interactions between lamina-other particles, and includes the nucleus deformation force when squeezing the model
 
-- utils: The folder contains two functions for calculating the contact prob between chromosome and nuclear landmarks (Lamina and Speckles), and the contact prob within chromosomes
+- utils: The folder contains functions for calculating the contact probabilities between chromosome and nuclear landmarks (Lamina and Speckles), and the contact probabilities among chromosomes. In addition, it also contains the coordinate transformation from OpenMM default unit to reduced unit, and the creation of the last configuration in one trajectory.

openNucleome/parameters/HFF_100KB/README.md

@@ -12,3 +12,5 @@ This folder includes all the parameter files needed for simulating a human nucle
 
 - types_param.txt: interaction parameters among various compartmental types
 
+- chromatin_info.txt: chromatin bead index (column 1), molecular info (column 2), and type info (column 3)
+

openNucleome/utils/coordinate_transformation.py

@@ -0,0 +1,29 @@
+import MDAnalysis as mda
+import numpy as np
+import sys
+
+def coor_transformation(traj_data, output_traj):
+    """
+    The function converts the trajectory with position in OpenMM default unit to reduced unit
+
+    Parameters
+    ----------
+    traj_data (string, required):
+        The trajectory file from a simulation
+
+    output_traj (string, required):
+        The trajectory with the reduced units
+    """
+
+    traj = mda.coordinates.LAMMPS.DCDReader(traj_data)
+
+    def scale_down_by_10(ts):
+        ts.positions /= 10
+        return ts
+
+    workflow = [scale_down_by_10]
+    traj.trajectory.add_transformations(*workflow)
+
+    with mda.Writer(output_traj, traj.n_atoms) as W:
+        for ts in traj.trajectory:
+            W.write(traj)

openNucleome/utils/final_frame.py

@@ -0,0 +1,64 @@
+import numpy as np
+import MDAnalysis as mda
+import sys
+import os
+
+current_path = os.path.abspath(__file__)
+father_path = os.path.abspath(os.path.dirname(current_path) + os.path.sep + '..')
+config_path = os.path.join(father_path, 'parameters', 'HFF_100KB')
+
+def final_frame(traj_data, type_final, output_config, n_chrom = 46, n_nuc = 300, n_spec = 1600, nuc_type = 5, lam_type = 8):
+    """
+    The function extracts the final frame of the previous simulation to restart a new simulation.
+
+    Parameters
+    ----------
+    traj_data (string, required):
+        The trajectory file from a simulation
+
+    type_final (string, required):
+        The file contains the type of each bead in the last frame
+
+    output_config (string, required):
+        The final frame of the trajectory
+
+    n_chrom (int, required):
+        The number of chromosomes (Default value = 46)
+
+    n_nuc (int, required):
+        The number of nucleolus (Default value = 300)
+
+    n_spec (int, required):
+        The number of speckles (Default value = 1600)
+
+    nuc_type (int, required):
+        The type index of nucleolus (Default value = 5)
+        
+    lam_type (int, required):
+        The type index of lamina (Default value = 8)
+    """
+
+    traj = mda.coordinates.LAMMPS.DCDReader(traj_data)
+    chr_info = np.loadtxt(os.path.join(config_path, 'chromatin_info.txt'))
+    X, zero = 'X', 0
+    type_final = np.loadtxt(type_final, dtype=int)
+    coor = traj.trajectory[-1].positions
+    n_all_beads = len(coor)
+    n_chrom_beads = len(chr_info)
+
+    file_write = open(output_config, 'w')
+    for i in range(n_chrom_beads):
+        if i == 0:
+            file_write.write('CRYST1   26.400   26.400   26.400  90.00  90.00  90.00 P 1           1\n')
+        file_write.write('ATOM%7s%3s%8s%4s%12.3f%8.3f%8.3f%6.2f%6.2f            \n'%(chr_info[i,0],chr_info[i,2],X,chr_info[i,1],coor[i,0],coor[i,1],coor[i,2],zero,zero))
+    for k in range(1,n_all_beads-n_chrom_beads+1):
+        idx_curr = n_chrom_beads+k
+        mol_curr = n_chrom+k
+        if k <= n_nuc:
+            file_write.write('ATOM%7s%3s%8s%4s%12.3f%8.3f%8.3f%6.2f%6.2f            \n'%(idx_curr,nuc_type,X,mol_curr,coor[idx_curr-1,0],coor[idx_curr-1,1],coor[idx_curr-1,2],zero,zero))
+        elif k <= n_nuc+n_spec:
+            file_write.write('ATOM%7s%3s%8s%4s%12.3f%8.3f%8.3f%6.2f%6.2f            \n'%(idx_curr,type_final[idx_curr-1],X,mol_curr,coor[idx_curr-1,0],coor[idx_curr-1,1],coor[idx_curr-1,2],zero,zero))
+        else:
+            file_write.write('ATOM%7s%3s%8s%4s%12.3f%8.3f%8.3f%6.2f%6.2f            \n'%(idx_curr,lam_type,X,mol_curr,coor[idx_curr-1,0],coor[idx_curr-1,1],coor[idx_curr-1,2],zero,zero))
+    file_write.write("END")
+    file_write.close()

tutorials/HFF_100KB/SphericalNucleus/BrownianDynamics/BrownianSimulation.ipynb

@@ -234,13 +234,13 @@
     "    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')"
+    "np.savetxt('results/compt_final_frame.txt', (np.array(model.compart_type)+1).reshape((-1,1)), fmt='%d')"
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -254,7 +254,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.8"
+   "version": "3.8.15"
   }
  },
  "nbformat": 4,