Implement initial support for bidirectional extrusion barriers

src/libmodle/cpu/scheduler_perturbate.cpp

@@ -469,7 +469,7 @@ void Simulation::simulate_window(Simulation::State& state, compressed_io::Writer
   for (usize i = 0; i < state.barriers.size(); ++i) {
     if (state.barriers.pos(i) >= state.deletion_begin &&
         state.barriers.pos(i) < state.deletion_begin + state.deletion_size) {
-      state.barriers.set(i, state.barriers.pos(i), state.barriers.direction(i), 0.0, 1.0,
+      state.barriers.set(i, state.barriers.pos(i), state.barriers.block_direction(i), 0.0, 1.0,
                          ExtrusionBarriers::State::INACTIVE);
       num_active_barriers--;
     }

src/libmodle/cpu/simulation.cpp

@@ -502,10 +502,10 @@ usize Simulation::release_lefs(const absl::Span<Lef> lefs, const ExtrusionBarrie
     assert(lefs[j].is_bound());
 
     auto is_lef_bar_hard_collision = [&](const auto c, dna::Direction d) {
-      assert(d == dna::REV || d == dna::FWD);
+      assert(d == dna::REV || d == dna::FWD || d == dna::BOTH);
       if (c.collision_occurred(CollisionT::LEF_BAR)) {
         assert(c.decode_index() <= barriers.size());
-        return barriers.direction(c.decode_index()) == d;
+        return barriers.block_direction(c.decode_index()) == d;
       }
       return false;
     };

src/libmodle/cpu/simulation_detect_collisions.cpp

@@ -182,8 +182,10 @@ void Simulation::detect_lef_bar_collisions(
 
     assert(j < lefs.size());
     // Probability of block is set based on the extr. barrier blocking direction
-    const auto& pblock = barriers.direction(i) == dna::REV ? this->lef_bar_major_collision_pblock
-                                                           : this->lef_bar_minor_collision_pblock;
+    const auto& pblock =
+        barriers.block_direction(i) == dna::REV || barriers.block_direction(i) == dna::BOTH
+            ? this->lef_bar_major_collision_pblock
+            : this->lef_bar_minor_collision_pblock;
 
     // Look for the first rev extr. unit that comes after the current barrier
     while (unit_pos <= barriers.pos(i)) {
@@ -229,8 +231,10 @@ void Simulation::detect_lef_bar_collisions(
     }
 
     // Probability of block is set based on the extr. barrier blocking direction
-    const auto& pblock = barriers.direction(i) == dna::FWD ? this->lef_bar_major_collision_pblock
-                                                           : this->lef_bar_minor_collision_pblock;
+    const auto& pblock =
+        barriers.block_direction(i) == dna::FWD || barriers.block_direction(i) == dna::BOTH
+            ? this->lef_bar_major_collision_pblock
+            : this->lef_bar_minor_collision_pblock;
     // Look for the next fwd unit that comes strictly before the current extr. barrier
     while (unit_pos >= barriers.pos(i)) {
       if (MODLE_UNLIKELY(--j == sentinel_idx)) {

src/libmodle/internal/extrusion_barriers.cpp

@@ -97,11 +97,11 @@ bp_t ExtrusionBarriers::pos(usize i) const noexcept {
 }
 const std::vector<bp_t>& ExtrusionBarriers::pos() const noexcept { return this->_pos; }
 
-dna::Direction ExtrusionBarriers::direction(usize i) const noexcept {
+dna::Direction ExtrusionBarriers::block_direction(usize i) const noexcept {
   this->assert_index_within_bounds(i);
   return this->_direction[i];
 }
-const std::vector<dna::Direction>& ExtrusionBarriers::direction() const noexcept {
+const std::vector<dna::Direction>& ExtrusionBarriers::block_direction() const noexcept {
   return this->_direction;
 }
 

src/libmodle/internal/extrusion_barriers_impl.hpp

@@ -13,6 +13,21 @@
 
 namespace modle {
 
+namespace internal {
+[[nodiscard]] constexpr dna::Direction char_to_strand(char d) {
+  switch (d) {
+    case '-':
+      return dna::REV;
+    case '+':
+      return dna::FWD;
+    case '.':
+      return dna::BOTH;
+    default:
+      throw std::runtime_error("Invalid DNA strand \"" + std::string{d} + "\"");
+  }
+}
+}  // namespace internal
+
 template <class BarrierIt, class StateIt>
 ExtrusionBarriers::ExtrusionBarriers(BarrierIt first_barrier, BarrierIt last_barrier,
                                      StateIt first_state, bool sort)
@@ -51,8 +66,9 @@ constexpr ExtrusionBarrier::ExtrusionBarrier(bp_t pos_, TP transition_prob_activ
     : pos(pos_),
       stp_active(transition_prob_active_to_active),
       stp_inactive(transition_prob_inactive_to_inactive),
-      blocking_direction(motif_direction == dna::FWD ? dna::REV : dna::FWD) {
-  assert(motif_direction == dna::FWD || motif_direction == dna::REV);
+      blocking_direction(motif_direction.complement()) {
+  assert(motif_direction == dna::FWD || motif_direction == dna::REV ||
+         motif_direction == dna::BOTH);
   assert(transition_prob_active_to_active() >= 0.0 && transition_prob_active_to_active() <= 1.0);
   assert(transition_prob_inactive_to_inactive() >= 0.0 &&
          transition_prob_inactive_to_inactive() <= 1.0);
@@ -64,8 +80,7 @@ constexpr ExtrusionBarrier::ExtrusionBarrier(bp_t pos_, TP transition_prob_activ
     : pos(pos_),
       stp_active(transition_prob_active_to_active),
       stp_inactive(transition_prob_inactive_to_inactive),
-      blocking_direction(motif_direction == '+' ? dna::REV : dna::FWD) {
-  assert(motif_direction == '+' || motif_direction == '-');
+      blocking_direction(internal::char_to_strand(motif_direction).complement()) {
   assert(transition_prob_active_to_active() >= 0.0 && transition_prob_active_to_active() <= 1.0);
   assert(transition_prob_inactive_to_inactive() >= 0.0 &&
          transition_prob_inactive_to_inactive() <= 1.0);
@@ -152,13 +167,24 @@ template <typename FormatContext>
 inline auto fmt::formatter<modle::ExtrusionBarrier>::format(const modle::ExtrusionBarrier& b,
                                                             FormatContext& ctx)
     -> decltype(ctx.out()) {
-  // ctx.out() is an output iterator to write to.
+  const auto strand = [&b]() {
+    if (b.blocking_direction == modle::dna::FWD) {
+      return "fwd";
+    }
+    if (b.blocking_direction == modle::dna::REV) {
+      return "rev";
+    }
+    if (b.blocking_direction == modle::dna::BOTH) {
+      return "both";
+    }
+  }();
+
   if (presentation == 's') {
-    return fmt::format_to(ctx.out(), "ExtrusionBarrier{{pos={}; motif_direction={}}}", b.pos,
-                          b.blocking_direction == modle::dna::FWD ? "rev" : "fwd");
+    return fmt::format_to(ctx.out(), "ExtrusionBarrier{{pos={}; blocking_direction={}}}", b.pos,
+                          strand);
   }
   assert(presentation == 'f');
-  return fmt::format_to(
-      ctx.out(), "ExtrusionBarrier{{pos={}; motif_direction={}; Pbb={:.4f}; Puu={:.4f}}}", b.pos,
-      b.blocking_direction == modle::dna::FWD ? "rev" : "fwd", b.stp_active, b.stp_inactive);
+  return fmt::format_to(ctx.out(),
+                        "ExtrusionBarrier{{pos={}; blocking_direction={}; Pbb={:.4f}; Puu={:.4f}}}",
+                        b.pos, strand, b.stp_active, b.stp_inactive);
 }

src/libmodle/internal/genome.cpp

@@ -409,11 +409,8 @@ usize Genome::import_barriers(absl::btree_set<Chromosome>& chromosomes,
                          "between 0 and 1, got {:.4g}."),
               record.chrom, record.chrom_start, record.chrom_end, record.score));
         }
-        if (record.strand == '+' || record.strand == '-') {
-          // For now we don't support barriers without strand information
-          chrom.add_extrusion_barrier(record, ctcf_prob_occ_to_occ, ctcf_prob_nocc_to_nocc);
-          ++tot_num_barriers;
-        }
+        chrom.add_extrusion_barrier(record, ctcf_prob_occ_to_occ, ctcf_prob_nocc_to_nocc);
+        ++tot_num_barriers;
       }
     }
     chrom.barriers().make_BST();

src/libmodle/internal/include/modle/extrusion_barriers.hpp

@@ -73,7 +73,7 @@ class ExtrusionBarriers {
   void clear() noexcept;
 
   [[nodiscard]] bp_t pos(usize i) const noexcept;
-  [[nodiscard]] dna::Direction direction(usize i) const noexcept;
+  [[nodiscard]] dna::Direction block_direction(usize i) const noexcept;
   [[nodiscard]] double stp_active(usize i) const noexcept;
   [[nodiscard]] double stp_inactive(usize i) const noexcept;
   [[nodiscard]] auto state(usize i) const noexcept -> State;
@@ -82,7 +82,7 @@ class ExtrusionBarriers {
   [[nodiscard]] double occupancy(usize i) const noexcept;
 
   [[nodiscard]] const std::vector<bp_t>& pos() const noexcept;
-  [[nodiscard]] const std::vector<dna::Direction>& direction() const noexcept;
+  [[nodiscard]] const std::vector<dna::Direction>& block_direction() const noexcept;
   [[nodiscard]] auto stp_active() const noexcept -> const std::vector<TP>&;
   [[nodiscard]] auto stp_inactive() const noexcept -> const std::vector<TP>&;
   [[nodiscard]] auto state() const noexcept -> const std::vector<State>&;
@@ -125,6 +125,11 @@ struct ExtrusionBarrier {
                                                             double occupancy) noexcept;
   static constexpr double compute_occupancy_from_stp(TP stp_active, TP stp_inactive) noexcept;
 };
+
+namespace internal {
+[[nodiscard]] constexpr dna::Direction char_to_strand(char d);
+}  // namespace internal
+
 }  // namespace modle
 
 template <>

src/modle/cli.cpp

@@ -118,7 +118,12 @@ static std::vector<CLI::App*> add_common_options(CLI::App& subcommand, modle::Co
       "interpreted as the extrusion barrier occupancy for the extrusion barrier described\n"
       "by the record.\n"
       "Barriers mapping on chromosomes not listed in the chrom.sizes file passed through\n"
-      "the --chrom-sizes option are ignored.")
+      "the --chrom-sizes option are ignored.\n"
+      "Currently, barriers are assumed to be either CTCF or generic bidirectional barriers.\n"
+      "If a barrier has + or - as its direction, that barrier is assumed to be a CTCF barrier,\n"
+      "and the provided strandness is interpreted as the motif/binding site direction.\n"
+      "Barriers with . as direction are treated as bidirectional barriers.\n"
+      "The current logic is a bit clunky. We have plans to improve this in the near future.")
       ->check(CLI::ExistingFile)
       ->required();
 

test/scripts/modle_integration_test_simple.sh

@@ -45,12 +45,24 @@ trap 'rm -rf -- "$outdir"' EXIT
                  -b <(xz -dc "$extr_barriers" | grep '^chr2') \
                  -o "$outdir/out" \
                  -r 20kb \
-                 --target-contact-density 20 \
-                 --ncells 4 \
-                 --max-burnin-epochs 5000
-
-# cp "$outdir/out.cool" /tmp/test/data/integration_tests/reference_001.cool
-echo "Comparing $outdir/out.cool with $data_dir/reference_001.cool..."
-h5diff --exclude-attribute / \
-       -c "$outdir/out.cool" \
-          "$data_dir/reference_001.cool"
+                 --contact-sampling-interval 5kb \
+                 --target-contact-density 5 \
+                 --ncells 2 \
+                 --max-burnin-epochs 750
+
+function compare_coolers() {
+  tgt="$1"
+  ref="$2"
+
+  # cp "$tgt" "$ref"
+  echo "Comparing $tgt with $ref..."
+  h5diff --exclude-attribute / \
+         -c "$tgt" "$ref"
+}
+
+if compare_coolers "$outdir/out.cool" "$data_dir/reference_001.cool"; then
+  >&2 echo "PASS!"
+else
+  >&2 echo "FAIL!"
+fi
+

test/units/simulation_cpu/simulation_simple_unit_test.cpp

@@ -587,7 +587,62 @@ TEST_CASE("Detect LEF-BAR collisions 002 - wo soft-collisions rev CTCFs",
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
-TEST_CASE("Detect LEF-BAR collisions 003 - w soft-collisions fwd CTCFs",
+TEST_CASE("Detect LEF-BAR collisions 003 - bidirectional barriers",
+          "[lef-bar-collisions][simulation][short]") {
+  const auto c = init_config(2, 2);
+  constexpr usize nlefs = 3;
+  constexpr usize nbarriers = 3;
+  auto rand_eng = DEFAULT_PRNG;
+
+  // clang-format off
+  std::array<Lef, nlefs> lefs{
+          construct_lef(0, 1, 0),
+          construct_lef(3, 4, 1),
+          construct_lef(5, 5, 2)
+  };
+  const auto barriers = construct_barriers(ExtrusionBarrier{2, 1.0, 0.0, '.'},
+                                           ExtrusionBarrier{4, 1.0, 0.0, '.'},
+                                           ExtrusionBarrier{8, 1.0, 0.0, '.'});
+  // clang-format on
+  REQUIRE(barriers.size() == nbarriers);
+
+  const std::array<usize, nlefs> rev_ranks{0, 1, 2};
+  const std::array<usize, nlefs> fwd_ranks{0, 1, 2};
+
+  std::array<bp_t, nlefs> rev_moves{0, 2, 2};
+  std::array<bp_t, nlefs> fwd_moves{2, 2, 2};
+
+  // clang-format off
+  const std::array<CollisionT, nlefs> rev_collisions_expected{CollisionT{},
+                                                              CollisionT{0, LEF_BAR},
+                                                              CollisionT{1, LEF_BAR}};
+  const std::array<CollisionT, nlefs> fwd_collisions_expected{CollisionT{0, LEF_BAR},
+                                                              CollisionT{},
+                                                              CollisionT{}};
+  // clang-format on
+
+  std::array<CollisionT, nlefs> rev_collisions;
+  std::array<CollisionT, nlefs> fwd_collisions;
+
+  const std::array<bp_t, nlefs> rev_moves_expected{0, 0, 0};
+  const std::array<bp_t, nlefs> fwd_moves_expected{0, 2, 2};
+
+  require_that_lefs_are_sorted_by_idx(lefs, rev_ranks, fwd_ranks);
+
+  Simulation{c, false}.test_detect_lef_bar_collisions(
+      lefs, rev_ranks, fwd_ranks, absl::MakeSpan(rev_moves), absl::MakeSpan(fwd_moves), barriers,
+      absl::MakeSpan(rev_collisions), absl::MakeSpan(fwd_collisions), rand_eng);
+  Simulation::test_correct_moves_for_lef_bar_collisions(lefs, barriers, absl::MakeSpan(rev_moves),
+                                                        absl::MakeSpan(fwd_moves), rev_collisions,
+                                                        fwd_collisions);
+
+  check_simulation_result(lefs, rev_moves, fwd_moves, rev_moves_expected, fwd_moves_expected,
+                          rev_collisions, rev_collisions_expected, fwd_collisions,
+                          fwd_collisions_expected);
+}
+
+// NOLINTNEXTLINE(readability-function-cognitive-complexity)
+TEST_CASE("Detect LEF-BAR collisions 004 - w soft-collisions fwd barriers",
           "[lef-bar-collisions][simulation][short]") {
   auto c = init_config(2, 2);
   c.lef_bar_major_collision_pblock = 1;
@@ -644,7 +699,7 @@ TEST_CASE("Detect LEF-BAR collisions 003 - w soft-collisions fwd CTCFs",
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
-TEST_CASE("Detect LEF-BAR collisions 004 - wo soft-collisions mixed CTCFs",
+TEST_CASE("Detect LEF-BAR collisions 005 - wo soft-collisions mixed barriers",
           "[lef-bar-collisions][simulation][short]") {
   const auto c = init_config(5, 5);
   const usize nlefs = 5;
@@ -707,8 +762,9 @@ TEST_CASE("Detect LEF-BAR collisions 004 - wo soft-collisions mixed CTCFs",
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
-TEST_CASE("Detect LEF-BAR collisions 005 - wo soft-collisions mixed CTCFs, different extr. speeds",
-          "[lef-bar-collisions][simulation][short]") {
+TEST_CASE(
+    "Detect LEF-BAR collisions 006 - wo soft-collisions mixed barriers, different extr. speeds",
+    "[lef-bar-collisions][simulation][short]") {
   const auto c = init_config(2, 5);
   const usize nlefs = 5;
   const usize nbarriers = 4;

test/units/simulation_internal/extrusion_barriers_test.cpp

@@ -32,34 +32,61 @@ TEST_CASE("Extrusion barriers - comparisons", "[barriers][simulation][short]") {
 TEST_CASE("Extrusion barriers - occupancy", "[barriers][simulation][short]") {
   const auto occupancy1 = 0.85;
   const auto occupancy2 = 0.93;
+  const auto occupancy3 = 0.35;
 
   const auto stp_inactive1 = 0.7;
   const auto stp_inactive2 = 0.65;
+  const auto stp_inactive3 = 0.9;
+
   const auto stp_active1 =
       ExtrusionBarrier::compute_stp_active_from_occupancy(stp_inactive1, occupancy1);
   const auto stp_active2 =
       ExtrusionBarrier::compute_stp_active_from_occupancy(stp_inactive2, occupancy2);
+  const auto stp_active3 =
+      ExtrusionBarrier::compute_stp_active_from_occupancy(stp_inactive3, occupancy3);
 
   using State = ExtrusionBarriers::State;
-  const auto barriers = ExtrusionBarriers{{100, 200},
-                                          {dna::FWD, dna::REV},
-                                          {stp_active1, stp_active2},
-                                          {stp_inactive1, stp_inactive2},
-                                          {State::ACTIVE, State::ACTIVE}};
+  const auto barriers = ExtrusionBarriers{{100, 200, 300},
+                                          {dna::FWD, dna::REV, dna::BOTH},
+                                          {stp_active1, stp_active2, stp_active3},
+                                          {stp_inactive1, stp_inactive2, stp_inactive3},
+                                          {State::ACTIVE, State::ACTIVE, State::ACTIVE}};
 
   CHECK(barriers.occupancy(0) == Catch::Approx(occupancy1));
   CHECK(barriers.occupancy(1) == Catch::Approx(occupancy2));
+  CHECK(barriers.occupancy(2) == Catch::Approx(occupancy3));
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Extrusion barriers - sort", "[barriers][simulation][short]") {
   using State = ExtrusionBarriers::State;
   auto barriers = ExtrusionBarriers{{100, 400, 200, 300},
-                                    {dna::FWD, dna::FWD, dna::FWD, dna::FWD},
-                                    {0.0, 0.0, 0.0, 0.0},
-                                    {0.0, 0.0, 0.0, 0.0},
+                                    {dna::FWD, dna::REV, dna::BOTH, dna::FWD},
+                                    {0.0, 0.1, 0.2, 0.3},
+                                    {0.3, 0.2, 0.1, 0.0},
                                     {State::ACTIVE, State::ACTIVE, State::ACTIVE, State::ACTIVE}};
 
+  const std::vector<bp_t> expected_pos{{100, 200, 300, 400}};
+  const std::vector<dna::Direction> expected_dir{{dna::FWD, dna::BOTH, dna::FWD, dna::REV}};
+  const std::vector<double> expected_stp_active{{0.0, 0.2, 0.3, 0.1}};
+  const std::vector<double> expected_stp_inactive{{0.3, 0.1, 0.0, 0.2}};
+
+  barriers.sort();
+  for (usize i = 0; i < barriers.size(); ++i) {
+    CHECK(barriers.pos(i) == expected_pos[i]);
+    CHECK(barriers.block_direction(i) == expected_dir[i]);
+    CHECK(barriers.stp_active(i) == expected_stp_active[i]);
+    CHECK(barriers.stp_inactive(i) == expected_stp_inactive[i]);
+  }
+
+  // Test sort w/ ties
+  barriers = ExtrusionBarriers{
+      {100, 400, 200, 300, 200, 100},
+      {dna::FWD, dna::FWD, dna::FWD, dna::FWD, dna::FWD, dna::FWD},
+      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+      {State::ACTIVE, State::ACTIVE, State::ACTIVE, State::ACTIVE, State::ACTIVE, State::ACTIVE}};
+
   barriers.sort();
   CHECK(std::is_sorted(barriers.pos().begin(), barriers.pos().end()));
 }