Linker length

README.md

@@ -3,32 +3,26 @@
 Physics-based simulator of chromatin architecture.
 
 ## Authors
-Joseph Wakim, Bruno Beltran, Andrew Spakowitz
+Joseph Wakim, Bruno Beltran, Angelika Hirsch, Andrew Spakowitz
 
 ## Quickstart
 
-Our physics-based chromatin simulator is publicly available in the [“Chromo” repository](https://github.com/SpakowitzLab/chromo) of the Spakowitz Lab GitHub account. Clone this repository to your local machine.
+The physics-based chromatin simulator is available in the [“Chromo” repository](https://github.com/AngelikaH123/chromo) of the Spakowitz Lab GitHub account. Clone this repository to your local machine.
 
-`git clone https://github.com/SpakowitzLab/chromo.git`
+`git clone git@github.com:AngelikaH123/chromo.git`
 
-Install the package dependencies of the simulator, which are listed in requirements.txt; we recommend that you do so in a separate virtual environment using Python 3.9.12. We demonstrate how to do so using the Conda package manager.
+Install the package dependencies and environment of the simulator by running the command:
 
-`conda create --name chromo python=3.9.12`
-
-`conda activate chromo`
-
-`pip install -r requirements.txt`
-
-Install our simulator package using pip by entering `pip install /path/to/root/directory` in the terminal, specifying the path to the root directory of the codebase on your local machine. If you would like to make changes to the codebase, use `pip install -e /path/to/root/directory` to install the simulator package in editable mode.
-
-`pip install /path/to/root/directory`
-
-`pip install -e /path/to/root/directory`
+'bash make_all.sh'
 
 During installation, all Cython code required by the Monte Carlo algorithm will be compiled. This may take several
-minutes. Once the package has been installed, the simulator is ready for use. Navigate to `simulations/examples` from the root directory to find example simulations.
+minutes. Once the package has been installed, the simulator is ready for use. Navigate to `chromo` from the root directory to find the run_simulation.py file for initializing the polymer object. Use the following command to run the file:
+
+'python run_simulation.py'
 
+This will generate the polymer object and compare it to a reference image, so the user can assess if the generated polymer is sufficiently similar to the control. 
 
+The design_document.docx file in the root directory provides additional information about the design of this project.  
 ## Abstract
 Chromatin organization is essential for differential gene expression in humans, enabling a diverse array of cell types to develop from the same genetic code. The architecture of chromatin is dictated by patterns of epigenetic marks–chemical modifications to DNA and histones–which must be reestablished during each cell replication cycle. Dysregulation of chromatin organization has drastic medical consequences and can contribute to aging, obesity, and cancer progression. During this study, we develop a Monte Carlo simulator called “Chromo” to investigate factors affecting chromatin organization. The simulator proposes and evaluates random configurational changes to chromatin, which include geometric transformations and effector protein binding/unbinding to the biopolymer. By iteratively sampling configurations, Chromo generates snapshots of energetically favorable chromatin architectures. Unlike previous models fit to experimental data, Chromo is rooted in fundamental polymer theory, allowing us to predict biophysical mechanisms governing chromatin organization. By leveraging a computationally efficient field theoretic approach, we can simulate full human chromosomes with nucleosome-scale resolution. To confirm the validity of our simulator, we reproduce theoretical chain statistics for flexible and semiflexible homopolymers and recapitulate heterochromatin compartments expected in chromatin contact maps. Using Chromo, we will evaluate how the quality of structural prediction changes with the simulation of additional epigenetic marks. We will also identify which combinations of marks best explain the architecture of the chromosome.
 

chromo/bead_position_analysis.py

@@ -0,0 +1,42 @@
+import numpy as np
+import sys
+import pandas as pd
+import matplotlib.pyplot as plt
+
+sherlock = False
+
+if sherlock:
+    file_path = "/scratch/users/ahirsch1/chromo_scratch/output/"
+    simulation_number = int(sys.argv[1])
+    snap_range1 = int(sys.argv[2])
+    snap_range1 = int(sys.argv[3])
+else:
+    file_path = "/Users/angelikahirsch/Documents/chromo/output/"
+    simulation_number = 218
+    snap_range1 = 177
+    snap_range2 = 179
+
+euclidean_distances = []
+for snap in range(snap_range1, snap_range2):
+    print(snap)
+    snap_path = f"{file_path}sim_{simulation_number}/poly_1-{snap}.csv"
+    df = pd.read_csv(snap_path)
+    coordinates = df[["('r', 'x')", "('r', 'y')", "('r', 'z')"]].values
+    num_beads = len(df)
+
+    for i in range(0, num_beads - 2):
+        difference = coordinates[i] - coordinates[i + 2]
+        distance = np.linalg.norm(difference)
+        euclidean_distances.append(distance/0.34)
+
+
+np.savetxt(f"{file_path}sim_{simulation_number}/{simulation_number}1,3bead_distances.txt", euclidean_distances)
+
+plt.hist(euclidean_distances, bins=20, color='blue', alpha=0.7)
+plt.xlabel('Distance in bp')
+plt.ylabel('Frequency')
+plt.title('Histogram of Distances between Beads Two Spaces Apart')
+plt.grid(True)
+plt.savefig(f"{file_path}sim_{simulation_number}/{simulation_number}1,3bead_distances.png")
+#plt.show()
+

chromo/beads.py

@@ -372,7 +372,7 @@ class Nucleosome(Bead):
 
     Attributes
     ----------
-    bead_length : double
+    bead_length : double ....changed to np.ndarray
         Spacing between the nucleosome and its neighbor
     rad : double
         Radius of spherical excluded volume around nucleosome
@@ -382,7 +382,7 @@ class Nucleosome(Bead):
 
     def __init__(
         self, id_: int, r: np.ndarray, *, t3: np.ndarray, t2: np.ndarray,
-        bead_length: float, rad: Optional[float] = 5,
+        bead_length: np.array, rad: Optional[float] = 5,
         states: Optional[np.ndarray] = None,
         binder_names: Optional[Sequence[str]] = None
     ):

chromo/mc/__init__.py

@@ -24,7 +24,8 @@
 from chromo.binders import ReaderProtein
 from chromo.binders import get_by_name, make_binder_collection
 from chromo.fields import UniformDensityField, FieldBase
-
+import chromo.twist_schedule as twist_schedule
+import chromo.util.temperature_schedule as temp_schedule
 
 F = TypeVar("F")    # Represents an arbitrary field
 STEPS = int         # Number of steps per MC save point
@@ -44,7 +45,9 @@ def _polymer_in_field(
     mc_move_controllers: Optional[List[Controller]] = None,
     random_seed: Optional[int] = 0,
     mu_schedule: Optional[Callable[[float], float]] = None,
-    output_dir: Optional[DIR] = '.',
+    lt_schedule: Optional[Callable[[str],float]] = None,
+    temperature_schedule: Optional[Callable[[str],float]] = None,
+    output_dir: Optional[DIR] = '.', # defaults to  current location, same level as chromo folder: name = output
     path_to_run_script: Optional[str] = None,
     path_to_chem_mods: Optional[List[str]] = None,
     run_command: Optional[str] = None,
@@ -126,17 +129,50 @@ def _polymer_in_field(
 
     t1_start = process_time()
     for mc_count in range(num_saves):
-
         # Simulated annealing
         if mu_schedule is not None:
             mu_adjust_factor = mu_schedule.function(mc_count, num_saves)
         else:
             mu_adjust_factor = 1
+        if temperature_schedule is not None:
+            if temperature_schedule == "linear decrease":
+                temperature_adjust_factor = temp_schedule.linear_decrease(mc_count, num_saves)
+            elif temperature_schedule == "logarithmic decrease":
+                temperature_adjust_factor = temp_schedule.logarithmic_decrease(mc_count, num_saves)
+            elif temperature_schedule == "decreasing stepwise":
+                temperature_adjust_factor = temp_schedule.decreasing_stepwise(mc_count, num_saves)
+            elif temperature_schedule == "sin decrease":
+                temperature_adjust_factor = temp_schedule.sin_decrease(mc_count, num_saves)
+            elif temperature_schedule == "no schedule":
+                temperature_adjust_factor = temp_schedule.no_schedule()
+            else:
+                print("Not a valid temperature schedule option")
+        else:
+            print("Missing lt schedule option")
+
+        if lt_schedule is not None:
+            if lt_schedule == "logarithmic increase":
+                lt_change = twist_schedule.logarithmic_increase(mc_count, num_saves)
+            elif lt_schedule == "exponential increase":
+                lt_change = twist_schedule.exponential_increase(mc_count, num_saves)
+            elif lt_schedule == "linear increase":
+                lt_change = twist_schedule.linear_increase(mc_count, num_saves)
+            elif lt_schedule == "increasing stepwise":
+                lt_change = twist_schedule.step_wise_increase(mc_count, num_saves)
+            elif lt_schedule == "increasing sawtooth":
+                lt_change = twist_schedule.increasing_sawtooth(mc_count, num_saves)
+            elif lt_schedule == "no schedule":
+                lt_change = polymers[0].lt  # set to the lt of the first polymer in the list
+            else:
+                print("Not a valid lt schedule option")
+        else:
+            print("Missing lt schedule option")
 
         decorator_timed_path(output_dir)(mc_sim)(
             polymers, binders, num_save_mc, mc_move_controllers, field,
-            mu_adjust_factor, random_seed
+            mu_adjust_factor, random_seed, temperature_adjust_factor, lt_change
         )
+
         for poly in polymers:
             poly.to_csv(
                 str(output_dir / Path(f"{poly.name}-{mc_count}.csv"))
@@ -149,6 +185,7 @@ def _polymer_in_field(
         print("Save point " + str(mc_count) + " completed")
 
     for polymer in polymers:
+        print("update log path " + str(output_dir))
         polymer.update_log_path(
             str(output_dir) + "/" + polymer.name + "_config_log.csv"
         )
@@ -196,6 +233,7 @@ def continue_polymer_in_field_simulation(
     random_seed : Optional[SEED]
         Random seed for replication of simulation (default = 0)
     """
+    # issue is somewhere here
     latest_output_subdir = get_latest_simulation(output_dir)
     latest_output_subdir_path = output_dir + "/" + latest_output_subdir
     polymer_names = find_polymers_in_output_dir(latest_output_subdir_path)

chromo/mc/mc_controller.py

@@ -254,8 +254,9 @@ def all_moves(
             move_amp_bounds=move_amp_bounds[move.__name__]
         ) for move in [
             mv.crank_shaft, mv.end_pivot, mv.slide, mv.tangent_rotation,
-            mv.change_binding_state
+            # mv.change_binding_state
         ]
+        #
     ]
     controllers[0].move.num_per_cycle = 30
     controllers[1].move.num_per_cycle = 1

chromo/mc/mc_sim.pxd

@@ -10,10 +10,11 @@ from chromo.fields cimport UniformDensityField as UDF
 cpdef void mc_sim(
     list polymers, readerproteins, long num_mc_steps,
     list mc_move_controllers, UDF field, double mu_adjust_factor,
-    long random_seed
+    long random_seed, double temperature_adjust_factor, double lt_value_adjust
 )
 
 cpdef void mc_step(
     MCAdapter adaptible_move, PolymerBase poly, readerproteins,
-    UDF field, bint active_field
-)
+    UDF field, bint active_field, double temperature_adjust_factor
+)
+

chromo/mc/mc_sim.pyx

@@ -25,12 +25,13 @@ from chromo.mc.moves import MCAdapter
 from chromo.mc.moves cimport MCAdapter
 from chromo.mc.mc_controller import Controller
 from chromo.fields cimport UniformDensityField as Udf
+from chromo.polymers import SSTWLC
 
 
 cpdef void mc_sim(
     list polymers, readerproteins, long num_mc_steps,
     list mc_move_controllers, Udf field, double mu_adjust_factor,
-    long random_seed
+    long random_seed, double temperature_adjust_factor, double lt_value_adjust
 ):
     """Perform Monte Carlo simulation.
 
@@ -75,32 +76,32 @@ cpdef void mc_sim(
     cdef long i, j, k, n_polymers
     cdef list active_fields
     cdef poly
-
     np.random.seed(random_seed)
-    n_polymers = len(polymers)
     if field.confine_type == "":
         active_fields = [poly.is_field_active() for poly in polymers]
     else:
         active_fields = [1 for _ in polymers]
-
     for k in range(num_mc_steps):
         for controller in mc_move_controllers:
             if controller.move.move_on == 1:
                 for j in range(controller.move.num_per_cycle):
                     for i in range(len(polymers)):
+                        polymers[i].lt = lt_value_adjust
+                        polymers[i]._find_parameters(polymers[i].bead_length) # coding best practice: accessing protected member of a class from outside the class?
                         poly = polymers[i]
-                        poly.mu_adjust_factor = mu_adjust_factor
                         active_field = active_fields[i]
+                        print("controller move")
+                        print(controller.move)
                         mc_step(
                             controller.move, poly, readerproteins, field,
-                            active_field
+                            active_field, temperature_adjust_factor
                         )
             controller.update_amplitudes()
 
 
 cpdef void mc_step(
     MCAdapter adaptible_move, PolymerBase poly, readerproteins,
-    Udf field, bint active_field
+    Udf field, bint active_field, double temperature_adjust_factor
 ):
     """Compute energy change and determine move acceptance.
 
@@ -126,16 +127,19 @@ cpdef void mc_step(
         Field affecting polymer in Monte Carlo step
     active_field: bool
         Indicator of whether or not the field is active for the polymer
+    temperature_adjust_factor: double
+        Divide the energy value used to select accepted moves by this factor
     """
     cdef double dE, exp_dE
     cdef int check_field = 0
     cdef long packet_size, n_inds
     cdef long[:] inds
 
+
+    #packet_size = 1 # will need to change back
     if poly in field and active_field:
         if adaptible_move.name != "tangent_rotation":
             check_field = 1
-
     packet_size = 20
     inds = adaptible_move.propose(poly)
     n_inds = len(inds)
@@ -161,14 +165,16 @@ cpdef void mc_step(
         elif dE < 0:
             exp_dE = 1
 
-    if (<double>rand() / RAND_MAX) < exp_dE:
+    if (<double>rand() / RAND_MAX) < exp_dE/temperature_adjust_factor: #temperature variation
+        #print(<double>rand() / RAND_MAX)
         adaptible_move.accept(
-            poly, dE, inds, n_inds, log_move=False, log_update=False
+            poly, dE, inds, n_inds, log_move=True, log_update=True
         )
         if check_field == 1:
             field.update_affected_densities()
 
     else:
         adaptible_move.reject(
-            poly, dE, log_move=False, log_update=False
+            poly, dE, log_move=True, log_update=True
         )
+

chromo/mc/moves.pyx

@@ -207,13 +207,15 @@ cdef class MCAdapter:
         self.acceptance_tracker.update_acceptance_rate(
             accept=1.0, log_update=log_update
         )
-        if log_move == 1:
+        if log_move == 1: # control whether acceptance trackers are on from here?
             self.acceptance_tracker.log_move(
                 self.amp_move,
                 self.amp_bead,
                 poly.last_amp_move,
                 poly.last_amp_bead,
-                dE
+                dE#,
+                #mc_stat.ConfigurationTracker(log_dir, sfe).RMSD
+               # mc_stat.ConfigurationTracker(log_path, )
             )
 
     cdef void reject(
@@ -242,7 +244,7 @@ cdef class MCAdapter:
         if log_move == 1:
             self.acceptance_tracker.log_move(
                 self.amp_move, self.amp_bead, poly.last_amp_move,
-                poly.last_amp_bead, dE
+                poly.last_amp_bead, dE#, mc_stat.ConfigurationTracker().RMSD
             )