Bed processing improvements

CMakeLists.txt

@@ -213,7 +213,7 @@ add_subdirectory(${CMAKE_CURRENT_BINARY_DIR}/googletest-src
                  ${CMAKE_CURRENT_BINARY_DIR}/googletest-build
                  EXCLUDE_FROM_ALL)
 
-
+# leviosam_test
 add_executable(leviosam_test src/leviosam_test.cpp)
 target_link_libraries(leviosam_test lvsam)
 target_link_libraries(leviosam_test ${HTS_LIB})
@@ -224,7 +224,7 @@ add_test(NAME leviosam_test COMMAND leviosam_test
     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/testdata
 )
 
-
+# chain_test
 add_executable(chain_test src/chain_test.cpp)
 target_link_libraries(chain_test lvsam)
 target_link_libraries(chain_test ${HTS_LIB})
@@ -235,6 +235,28 @@ add_test(NAME chain_test COMMAND chain_test
     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/testdata
 )
 
+# bed_test
+add_executable(bed_test src/bed_test.cpp)
+target_link_libraries(bed_test lvsam)
+target_link_libraries(bed_test ${HTS_LIB})
+target_link_libraries(bed_test ${SDSL_LIB})
+target_link_libraries(bed_test gtest gtest_main)
+
+add_test(NAME bed_test COMMAND bed_test
+    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/testdata
+)
+
+# reconcile_test
+add_executable(reconcile_test src/reconcile_test.cpp)
+target_link_libraries(reconcile_test lvsam)
+target_link_libraries(reconcile_test ${HTS_LIB})
+target_link_libraries(reconcile_test ${SDSL_LIB})
+target_link_libraries(reconcile_test gtest gtest_main)
+
+add_test(NAME reconcile_test COMMAND reconcile_test
+    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/testdata
+)
+
 INSTALL(TARGETS leviosam2 DESTINATION bin)
 INSTALL(FILES src/leviosam.hpp DESTINATION include)
 

src/IITree.h

@@ -89,15 +89,22 @@ class IITree {
                 size_t i, i0 = z.x >> z.k << z.k, i1 = i0 + (1LL<<(z.k+1)) - 1;
                 if (i1 >= a.size()) i1 = a.size();
                 for (i = i0; i < i1 && a[i].st < en; ++i)
-                    if (st < a[i].en) // if overlap, append to out[]
-                        out.push_back(i);
+                    if (st < a[i].en) { // if overlap, append to out[]
+                        int _en = std::min(en, a[i].en);
+                        int _st = std::max(st, a[i].st);
+                        out.push_back(std::max(1, _en - _st));
+                    }
             } else if (z.w == 0) { // if left child not processed
                 size_t y = z.x - (1LL<<(z.k-1)); // the left child of z.x; NB: y may be out of range (i.e. y>=a.size())
                 stack[t++] = StackCell(z.k, z.x, 1); // re-add node z.x, but mark the left child having been processed
                 if (y >= a.size() || a[y].max > st) // push the left child if y is out of range or may overlap with the query
                     stack[t++] = StackCell(z.k - 1, y, 0);
             } else if (z.x < a.size() && a[z.x].st < en) { // need to push the right child
-                if (st < a[z.x].en) out.push_back(z.x); // test if z.x overlaps the query; if yes, append to out[]
+                if (st < a[z.x].en) {
+                    int _en = std::min(en, a[z.x].en);
+                    int _st = std::max(st, a[z.x].st);
+                    out.push_back(std::max(1, _en - _st)); // test if z.x overlaps the query; if yes, append to out[]
+                }
                 stack[t++] = StackCell(z.k - 1, z.x + (1LL<<(z.k-1)), 0); // push the right child
             }
         }

src/bed.cpp

@@ -33,7 +33,6 @@ void Bed::init(const std::string &fn) {
         std::string line;
         while (std::getline(bed_f, line)) {
             cnt += 1;
-            // std::cerr << line << "\n";
             if (add_interval(line) == false) {
                 std::cerr << "[E::Bed] Failed to update BED record " << line
                           << "\n";
@@ -46,9 +45,10 @@ void Bed::init(const std::string &fn) {
               << " contigs) from " << fn << "\n";
 }
 
+/* Index the interval trees for each contig and returns the number of contigs
+ */
 int Bed::index() {
     int contig_cnt = 0;
-    // Index the interval trees for each contig
     for (auto &itr : intervals) {
         itr.second.index();
         contig_cnt += 1;
@@ -76,14 +76,36 @@ bool Bed::add_interval(const std::string &line) {
     return true;
 }
 
-// Interval intersect query
+/* Interval intersect query
+ * Args:
+ *   isec_threshold:
+ *    - <0: return true if overlap >= 1bp
+ *    - 1>= `isec_threshold` > 0: return true if overlap fraction >= `isec_threshold`
+ *    - `isec_threshold` > 1: return true if overlap size >= `isec_threshold` (bp)
+*/
 bool Bed::intersect(const std::string &contig, const size_t &pos1,
-                    const size_t &pos2) {
+                    const size_t &pos2, const float &isec_threshold) {
     if (!is_valid) return false;
+    if (pos2 <= pos1) return false;
     if (intervals.find(contig) == intervals.end()) return false;
     std::vector<size_t> isec;
-    intervals[contig].overlap(pos1, pos2, isec);
-    return (isec.size() > 0);
+    bool ovlp_status = intervals[contig].overlap(pos1, pos2, isec);
+    if (!ovlp_status) return false;
+    // <=0 means any overlap is good
+    if (isec_threshold <= 0) return (isec.size() > 0);
+
+    int max_overlap = 0;
+    for (auto& _o: isec) {
+        if (_o > max_overlap) max_overlap = _o;
+    }
+    if (isec_threshold > 0 && isec_threshold <= 1)
+        return (max_overlap >= (pos2 - pos1) * isec_threshold);
+
+    if (isec_threshold > 1) 
+        return (max_overlap >= isec_threshold);
+
+    return false;
+    // return (isec.size() > 0);
 }
 
 // Point intersect query

src/bed.hpp

@@ -32,7 +32,7 @@ class Bed {
     int index();
     bool add_interval(const std::string &line);
     bool intersect(const std::string &contig, const size_t &pos1,
-                   const size_t &pos2);
+                   const size_t &pos2, const float &isec_threshold = 0);
     bool intersect(const std::string &contig, const size_t &pos);
     BedMap get_intervals();
     std::string get_fn();

src/bed_test.cpp

@@ -0,0 +1,89 @@
+/*
+ * bed_test.cpp
+ *
+ * Test bed related functions for levioSAM2
+ *
+ * Author: Nae-Chyun Chen
+ * Dept. of Computer Science, Johns Hopkins University
+ *
+ * Distributed under the MIT license
+ * https://github.com/milkschen/leviosam2
+ */
+#include "bed.hpp"
+
+#include <unistd.h>
+
+#include <iostream>
+
+#include "gtest/gtest.h"
+#include "leviosam_utils.hpp"
+
+/* Bed tests */
+TEST(BedTest, Simple) {
+    std::string bed_fn("test_small.bed");
+    BedUtils::Bed bed(bed_fn);
+    EXPECT_EQ(bed.get_fn(), bed_fn);
+    EXPECT_EQ(bed.index(), 3);
+    EXPECT_EQ(bed.intersect("chr1", 450, 499), false);
+    EXPECT_EQ(bed.intersect("chr1", 450, 501), true);
+}
+
+TEST(BedTest, AddInterval) {
+    std::string bed_fn("test_small.bed");
+    BedUtils::Bed bed(bed_fn);
+    EXPECT_EQ(bed.intersect("chr1", 55), false);
+    EXPECT_EQ(bed.index(), 3);
+    EXPECT_EQ(bed.add_interval("chr1	50	60"), true);
+    // chr1 already exists, so cnt is still 3
+    EXPECT_EQ(bed.index(), 3);
+    EXPECT_EQ(bed.intersect("chr1", 55), true);
+    EXPECT_EQ(bed.intersect("chr1", 65), false);
+    EXPECT_EQ(bed.intersect("chr1", 55, 58), true);
+    EXPECT_EQ(bed.intersect("chr1", 45, 80), true);
+    EXPECT_EQ(bed.intersect("chr1", 45, 55), true);
+}
+
+TEST(BedTest, IntersectWithFrac) {
+    std::string bed_fn("test_small.bed");
+    BedUtils::Bed bed(bed_fn);
+    EXPECT_EQ(bed.intersect("chr1", 450, 650, 1.0), false);
+    EXPECT_EQ(bed.intersect("chr1", 550, 580, 1.0), true);
+    EXPECT_EQ(bed.intersect("chr1", 550, 580, 0.99), true);
+    EXPECT_EQ(bed.intersect("chr1", 550, 700, 0.3), true);
+    EXPECT_EQ(bed.intersect("chr1", 550, 700, 0.35), false);
+    EXPECT_EQ(bed.intersect("chr1", 550, 700, 30), true);
+    EXPECT_EQ(bed.intersect("chr1", 550, 580, 0), true);
+    EXPECT_EQ(bed.intersect("chr1", 550, 580, -1), true);
+}
+
+TEST(BedTest, AddBedRecord) {
+    auto b = BedUtils::Bed();
+    b.add_interval("chr1	972	1052");
+    b.add_interval("chr1	1053	1085");
+    b.add_interval("chr1	1150	1198");
+    b.add_interval("chr1	1199	1262");
+    b.add_interval("chr2	340	498");
+    b.add_interval("chr2	599	1000");
+    int idx_contig_cnt = b.index();
+    int contig_cnt = 0;
+    // Index the interval trees for each contig
+    for (auto& itr: b.get_intervals()) {
+        contig_cnt += 1;
+    }
+    EXPECT_EQ(contig_cnt, 2);
+    EXPECT_EQ(contig_cnt, idx_contig_cnt);
+    EXPECT_EQ(b.intersect("chr1", 1000), 1);
+    EXPECT_EQ(b.intersect("chr1", 1052), 0);
+    EXPECT_EQ(b.intersect("chr2", 339), 0);
+    EXPECT_EQ(b.intersect("chr2", 340), 1);
+    EXPECT_EQ(b.intersect("chr2", 497), 1);
+    EXPECT_EQ(b.intersect("chr2", 498), 0);
+    EXPECT_EQ(b.intersect("chr2", 340), 1);
+    EXPECT_EQ(b.intersect("chr2", 340), 1);
+    EXPECT_EQ(b.intersect("chr2", 598), 0);
+    EXPECT_EQ(b.intersect("chr2", 599), 1);
+    EXPECT_EQ(b.intersect("chr2", 600), 1);
+    EXPECT_EQ(b.intersect("chr1", 1500), 0);
+    EXPECT_EQ(b.intersect("chr15", 1000), 0);
+}
+