Merge pull request #268 from CovertLab/superhelical-fix

Remove molecules at same coordinates as replisomes

doc/ci.rst

@@ -37,10 +37,19 @@ tasks are:
 Logs from these tests can be viewed on the GitHub website and we strongly
 recommend that you get all of these tests passing before merging a PR.
 
-Additionally, we would appreciate if you added tests to improve our test coverage
-and improve the likelihood that we catch bugs. This could be as simple as a Python
-function with the ``test_`` prefix that ensures some bit of code changed in your
-PR works as intended with a few test cases.
+When you submit a pull request (PR), a bot will comment with a table showing the current code
+coverage of the ``pytest`` suite. As of December 2024, coverage is less than 30%
+and even that is misleading because many of the tests are "migration" tests
+that compare vEcoli to a snapshot of the original whole cell model. These tests will
+be removed in the future as vEcoli is further developed. Additionally, these tests do
+not tell us whether the code is working as intended, only that it is working the same
+way it worked in the original model. Ideally, we would like to increase test coverage
+by adding unit tests which actually test edge cases and ensure the code does what it
+is supposed to do.
+
+To that end, we would appreciate if you added tests as part of your pull requests.
+This could be as simple as a Python function with the ``test_`` prefix that ensures
+the code added or modified in your PR works as intended using a few test cases.
 
 -------
 Jenkins

doc/workflows.rst

@@ -740,9 +740,8 @@ the output directory specified via ``out_dir`` or ``out_uri`` under the
         - ``nextflow_workdirs``: Contains all working directories for Nextflow jobs.
           Required for resume functionality described in :ref:`fault_tolerance`. Can
           also go to work directory for a job (consult files described in :ref:`progress`
-          or ``{experiment ID}_report.html``) and run ``bash .command.sh`` with
-          breakpoints set in the relevant code (``import ipdb; ipdb.set_trace()``)
-          for debugging.
+          or ``{experiment ID}_report.html``) for debugging. See :ref:`make_and_test`
+          for more information.
 
 .. tip::
   To save space, you can safely delete ``nextflow_workdirs`` after you are finished
@@ -795,6 +794,8 @@ in a workflow called ``agitated_mendel``::
   nextflow log agitated_mendel -f name,stderr,workdir -F "status == 'FAILED'"
 
 
+.. _make_and_test:
+
 Make and Test Fixes
 ===================
 
@@ -810,11 +811,14 @@ Add breakpoints to any Python file with the following line::
   import ipdb; ipdb.set_trace()
 
 Then, navigate to the working directory (see :ref:`troubleshooting`) for a
-failing process. ``bash .command.run`` should re-run the job and pause upon
-reaching the breakpoints you set. You should now be in an ipdb shell which
-you can use to examine variable values or step through the code.
+failing process. Invoke ``uv run --env-file {} .command.run``, replacing
+the curly braces with the path to the ``.env`` file in your cloned repository.
+This should re-run the job and pause upon reaching the breakpoints you set.
+You should now be in an ipdb shell which you can use to examine variable values
+or step through the code.
 
 After fixing the issue, you can resume the workflow (avoid re-running
 already successful jobs) by navigating back to the directory in which you
-originally started the workflow and issuing the same command with the
-``--resume`` option (see :ref:`fault_tolerance`).
+originally started the workflow and issuing the same command
+(:py:mod:`runscripts.workflow`) with the ``--resume`` option
+(see :ref:`fault_tolerance`).

ecoli/composites/ecoli_master_tests.py

@@ -149,8 +149,8 @@ def test_division_topology():
         "ignore",
         message="Incompatible schema "
         "assignment at .+ Trying to assign the value <bound method "
-        "UniqueNumpyUpdater\.updater .+ to key updater, which already "
-        "has the value <bound method UniqueNumpyUpdater\.updater",
+        r"UniqueNumpyUpdater\.updater .+ to key updater, which already "
+        r"has the value <bound method UniqueNumpyUpdater\.updater",
     )
     sim.ecoli_experiment = Engine(**experiment_config)
 

ecoli/experiments/ecoli_engine_process.py

@@ -460,8 +460,8 @@ def run_simulation(config):
         "ignore",
         message="Incompatible schema "
         "assignment at .+ Trying to assign the value <bound method "
-        "UniqueNumpyUpdater\.updater .+ to key updater, which already "
-        "has the value <bound method UniqueNumpyUpdater\.updater",
+        r"UniqueNumpyUpdater\.updater .+ to key updater, which already "
+        r"has the value <bound method UniqueNumpyUpdater\.updater",
     )
     engine = Engine(
         processes=composite.processes,

ecoli/experiments/ecoli_master_sim.py

@@ -820,8 +820,8 @@ def run(self):
             "ignore",
             message="Incompatible schema "
             "assignment at .+ Trying to assign the value <bound method "
-            "UniqueNumpyUpdater.updater .+ to key updater, which already "
-            "has the value <bound method UniqueNumpyUpdater.updater",
+            r"UniqueNumpyUpdater\.updater .+ to key updater, which already "
+            r"has the value <bound method UniqueNumpyUpdater\.updater",
         )
         self.ecoli_experiment = Engine(**experiment_config)
 

ecoli/library/schema.py

@@ -126,7 +126,7 @@ def __array_finalize__(self, obj):
         # Views should inherit metadata from parent
         self.metadata = getattr(obj, "metadata", None)
 
-    def __array_wrap__(self, out_arr, context=None):
+    def __array_wrap__(self, out_arr, context=None, return_scalar=False):
         # If the result is a scalar, return it as a base scalar type
         if out_arr.shape == ():
             return out_arr.item()

ecoli/processes/chromosome_structure.py

@@ -12,14 +12,20 @@
 import numpy.typing as npt
 import warnings
 from vivarium.core.process import Step
+from vivarium.core.composer import Composer
+from vivarium.core.engine import Engine
 
+from ecoli.processes.global_clock import GlobalClock
+from ecoli.processes.unique_update import UniqueUpdate
 from ecoli.processes.registries import topology_registry
 from ecoli.library.schema import (
     listener_schema,
     numpy_schema,
     attrs,
     bulk_name_to_idx,
+    get_free_indices,
 )
+from ecoli.library.json_state import get_state_from_file
 from wholecell.utils.polymerize import buildSequences
 
 # Register default topology for this process, associating it with process name
@@ -64,15 +70,15 @@ class ChromosomeStructure(Step):
     topology = TOPOLOGY
     defaults = {
         # Load parameters
-        "RNA_sequences": [],
+        "rna_sequences": [],
         "protein_sequences": [],
         "n_TUs": 1,
         "n_TFs": 1,
         "n_amino_acids": 1,
         "n_fragment_bases": 1,
         "replichore_lengths": [0, 0],
         "relaxed_DNA_base_pairs_per_turn": 1,
-        "terC_index": 1,
+        "terC_index": -1,
         "calculate_superhelical_densities": False,
         # Get placeholder value for chromosome domains without children
         "no_child_place_holder": -1,
@@ -87,6 +93,13 @@ class ChromosomeStructure(Step):
         "water": "water",
         "seed": 0,
         "emit_unique": False,
+        "rna_ids": [],
+        "n_mature_rnas": 0,
+        "mature_rna_ids": [],
+        "mature_rna_end_positions": [],
+        "mature_rna_nt_counts": [],
+        "unprocessed_rna_index_mapping": {},
+        "time_step": 1.0,
     }
 
     # Constructor
@@ -307,11 +320,15 @@ def get_removed_molecules_mask(domain_indexes, coordinates):
                     ]
 
                     # Get mask for molecules on this domain that are out of range
+                    # It's rare but we have to remove molecules at the exact same
+                    # coordinates as the replisomes as well so that they do not break
+                    # the chromosome segment calculations if they are removed by a
+                    # different process (hence, >= and <= instead of > and <)
                     domain_mask = np.logical_and.reduce(
                         (
                             domain_indexes == domain_index,
-                            coordinates > domain_replisome_coordinates.min(),
-                            coordinates < domain_replisome_coordinates.max(),
+                            coordinates >= domain_replisome_coordinates.min(),
+                            coordinates <= domain_replisome_coordinates.max(),
                         )
                     )
 
@@ -1109,3 +1126,154 @@ def _compute_new_segment_attributes(
             "boundary_coordinates": new_boundary_coordinates,
             "linking_numbers": new_linking_numbers,
         }
+
+
+def test_superhelical_bug():
+    """
+    Test that ChromosomeStructure correctly handles edge case where RNAP and replisome
+    share the same genomic coordinates.
+    """
+    # Get topology for UniqueUpdate Steps
+    unique_topology = TOPOLOGY.copy()
+    for non_unique in [
+        "bulk",
+        "listeners",
+        "global_time",
+        "timestep",
+        "next_update_time",
+    ]:
+        unique_topology.pop(non_unique)
+
+    class TestComposer(Composer):
+        def generate_processes(self, config):
+            return {
+                "chromosome_structure": ChromosomeStructure(
+                    {
+                        "inactive_RNAPs": "APORNAP-CPLX[c]",
+                        "ppi": "PPI[c]",
+                        "active_tfs": ["CPLX-125[c]"],
+                        "ribosome_30S_subunit": "CPLX0-3953[c]",
+                        "ribosome_50S_subunit": "CPLX0-3962[c]",
+                        "amino_acids": ["L-ALPHA-ALANINE[c]"],
+                        "water": "WATER[c]",
+                        "mature_rna_ids": ["alaT-tRNA[c]"],
+                        "fragmentBases": ["polymerized_ATP[c]"],
+                        "replichore_lengths": [2000, 2000],
+                        "calculate_superhelical_densities": True,
+                    }
+                ),
+                "unique_update": UniqueUpdate({"unique_topo": unique_topology}),
+                "global_clock": GlobalClock(),
+            }
+
+        def generate_topology(self, config):
+            return {
+                "rnap_remover": {
+                    "active_RNAPs": ("unique", "active_RNAP"),
+                    "global_time": ("global_time",),
+                },
+                "chromosome_structure": TOPOLOGY,
+                "unique_update": unique_topology,
+                "global_clock": {
+                    "global_time": ("global_time",),
+                    "next_update_time": ("next_update_time",),
+                },
+            }
+
+        def generate_flow(self, config):
+            return {
+                "chromosome_structure": [],
+                "unique_update": [("chromosome_structure",)],
+            }
+
+    composer = TestComposer()
+    template_initial_state = get_state_from_file()
+    # Zero out all unique molecules
+    for unique_mol in template_initial_state["unique"].values():
+        unique_mol.flags.writeable = True
+        unique_mol["_entryState"] = 0
+        unique_mol.flags.writeable = False
+    # Set up chromosome domain 1
+    full_chromosomes = template_initial_state["unique"]["full_chromosome"]
+    full_chromosomes.flags.writeable = True
+    full_chromosomes["_entryState"][0] = 1
+    full_chromosomes["domain_index"][0] = 1
+    full_chromosomes.flags.writeable = False
+    chromosome_domains = template_initial_state["unique"]["chromosome_domain"]
+    chromosome_domains.flags.writeable = True
+    chromosome_domains["_entryState"][0] = 1
+    chromosome_domains["domain_index"][0] = 1
+    chromosome_domains["child_domains"][0] = [-1, -1]
+    chromosome_domains.flags.writeable = False
+    # Add two replisomes at -1000 and 1000
+    active_replisomes = template_initial_state["unique"]["active_replisome"]
+    active_replisomes.flags.writeable = True
+    active_replisomes["_entryState"][:2] = 1
+    active_replisomes["domain_index"][:2] = 1
+    active_replisomes["coordinates"][:2] = [-1000, 1000]
+    active_replisomes["unique_index"][:2] = [1, 2]
+    active_replisomes.flags.writeable = False
+    # Add four RNAPs: two that coincide with replisomes, two that do not
+    active_RNAP = template_initial_state["unique"]["active_RNAP"]
+    active_RNAP.flags.writeable = True
+    active_RNAP["_entryState"][:4] = 1
+    active_RNAP["domain_index"][:4] = 1
+    active_RNAP["coordinates"][:4] = [-1500, -1000, 1000, 1500]
+    active_RNAP["is_forward"][:4] = [True, True, True, True]
+    active_RNAP["unique_index"][:4] = [3, 4, 5, 6]
+    active_RNAP.flags.writeable = False
+    # Add chromosomal segments between replisome and RNAPs on either side
+    # of origin of replication (coordinates offset from current timestep)
+    chromosomal_segments, _ = get_free_indices(
+        template_initial_state["unique"]["chromosomal_segment"], 2
+    )
+    chromosomal_segments.flags.writeable = True
+    chromosomal_segments["_entryState"][:2] = 1
+    chromosomal_segments["boundary_coordinates"][:2] = [[-1400, -800], [800, 1400]]
+    chromosomal_segments["boundary_molecule_indexes"][:2] = [[3, 1], [2, 6]]
+    chromosomal_segments["domain_index"][:2] = 1
+    chromosomal_segments["linking_number"][:2] = 1
+    chromosomal_segments.flags.writeable = False
+    template_initial_state["unique"]["chromosomal_segment"] = chromosomal_segments
+    # Since unique numpy updater is an class method, internal
+    # deepcopying in vivarium-core causes this warning to appear
+    warnings.filterwarnings(
+        "ignore",
+        message="Incompatible schema "
+        "assignment at .+ Trying to assign the value <bound method "
+        r"UniqueNumpyUpdater\.updater .+ to key updater, which already "
+        r"has the value <bound method UniqueNumpyUpdater\.updater",
+    )
+    engine = Engine(
+        composite=composer.generate(),
+        initial_state=template_initial_state,
+    )
+    engine.update(1)
+    state = engine.state.get_value()
+    # Check that conflicting RNAPs were removed
+    active_RNAPs = state["unique"]["active_RNAP"][
+        state["unique"]["active_RNAP"]["_entryState"].view(np.bool_)
+    ]
+    assert len(active_RNAPs) == 2
+    assert np.all(active_RNAPs["domain_index"] == 1)
+    assert np.all(active_RNAPs["coordinates"] == np.array([-1500, 1500]))
+    assert np.all(active_RNAPs["unique_index"] == np.array([3, 6]))
+    # Check that chromosomal segments were added for the span between
+    # each replisome and terC in addition to the existing ones between
+    # replisomes and non-conflicting RNAPs
+    active_chromosome_segments = state["unique"]["chromosomal_segment"][
+        state["unique"]["chromosomal_segment"]["_entryState"].view(np.bool_)
+    ]
+    assert len(active_chromosome_segments) == 4
+    assert np.all(
+        active_chromosome_segments["boundary_coordinates"]
+        == np.array([[-2000, -1500], [-1500, -1000], [1000, 1500], [1500, 2000]])
+    )
+    assert np.all(
+        active_chromosome_segments["boundary_molecule_indexes"]
+        == np.array([[-1, 3], [3, 1], [2, 6], [6, -1]])
+    )
+
+
+if __name__ == "__main__":
+    test_superhelical_bug()

ecoli/processes/engine_process.py

@@ -289,8 +289,8 @@ def __init__(self, parameters=None):
             "ignore",
             message="Incompatible schema "
             "assignment at .+ Trying to assign the value <bound method "
-            "UniqueNumpyUpdater\.updater .+ to key updater, which already "
-            "has the value <bound method UniqueNumpyUpdater\.updater",
+            r"UniqueNumpyUpdater\.updater .+ to key updater, which already "
+            r"has the value <bound method UniqueNumpyUpdater\.updater",
         )
         self.sim = Engine(
             processes=processes,

ecoli/processes/rna_interference.py

@@ -213,8 +213,8 @@ def test_rna_interference(return_data=False):
         "ignore",
         message="Incompatible schema "
         "assignment at .+ Trying to assign the value <bound method "
-        "UniqueNumpyUpdater.updater .+ to key updater, which already "
-        "has the value <bound method UniqueNumpyUpdater.updater",
+        r"UniqueNumpyUpdater\.updater .+ to key updater, which already "
+        r"has the value <bound method UniqueNumpyUpdater\.updater",
     )
     experiment = Engine(
         processes={"rna-interference": rna_inter},

ecoli/processes/transcript_elongation.py

@@ -751,8 +751,8 @@ def make_elongation_rates(random, base, time_step, variable_elongation=False):
         "ignore",
         message="Incompatible schema "
         "assignment at .+ Trying to assign the value <bound method "
-        "UniqueNumpyUpdater.updater .+ to key updater, which already "
-        "has the value <bound method UniqueNumpyUpdater.updater",
+        r"UniqueNumpyUpdater\.updater .+ to key updater, which already "
+        r"has the value <bound method UniqueNumpyUpdater\.updater",
     )
     engine = Engine(**settings, initial_state=deepcopy(initial_state))
     engine.run_for(100)
@@ -813,8 +813,8 @@ def make_elongation_rates(random, base, time_step, variable_elongation=False):
         "ignore",
         message="Incompatible schema "
         "assignment at .+ Trying to assign the value <bound method "
-        "UniqueNumpyUpdater.updater .+ to key updater, which already "
-        "has the value <bound method UniqueNumpyUpdater.updater",
+        r"UniqueNumpyUpdater\.updater .+ to key updater, which already "
+        r"has the value <bound method UniqueNumpyUpdater\.updater",
     )
     engine = Engine(**settings, initial_state=deepcopy(initial_state))
     engine.run_for(100)

ecoli/processes/unique_update.py

@@ -3,8 +3,8 @@
 
 
 class UniqueUpdate(Step):
-    """Placed after before and after all Steps to ensure that
-    unique molecules are completely up-to-date"""
+    """Placed after all Steps of each execution layer (see :ref:`partitioning`)
+    to ensure that unique molecules are completely up-to-date"""
 
     name = "unique-update"
 

reconstruction/ecoli/dataclasses/process/rna_decay.py

@@ -8,6 +8,7 @@
 from autograd import jacobian
 import numpy as np
 
+
 class RnaDecay(object):
     """RnaDecay"""