BED feature

.gitmodules

@@ -19,4 +19,7 @@
 [submodule "vendor/GenomeWorks"]
 	path = vendor/GenomeWorks
 	url = https://github.com/clara-parabricks/GenomeWorks.git
+[submodule "vendor/intervaltree"]
+	path = vendor/intervaltree
+	url = https://github.com/ekg/intervaltree.git
 	branch = master

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.2)
 project(racon)
-set(racon_version 1.4.17)
+set(racon_version 1.5.0)
 
 set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib)
 set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib)
@@ -27,13 +27,16 @@ if(racon_enable_cuda)
 endif()
 
 include_directories(${PROJECT_SOURCE_DIR}/src)
+include_directories(${PROJECT_SOURCE_DIR}/vendor/intervaltree)
 
 set(racon_sources
     src/main.cpp
+    src/bed.cpp
     src/logger.cpp
     src/polisher.cpp
     src/overlap.cpp
     src/sequence.cpp
+    src/util.cpp
     src/window.cpp)
 
 if(racon_enable_cuda)
@@ -103,11 +106,15 @@ if (racon_build_tests)
     include_directories(${PROJECT_SOURCE_DIR}/src)
 
     set(racon_test_sources
+        test/bed_test.cpp
         test/racon_test.cpp
+        test/util_test.cpp
+        src/bed.cpp
         src/logger.cpp
         src/polisher.cpp
         src/overlap.cpp
         src/sequence.cpp
+        src/util.cpp
         src/window.cpp)
 
     if (racon_enable_cuda)

meson.build

@@ -1,7 +1,7 @@
 project(
   'Racon',
   'cpp',
-  version : '1.4.13',
+  version : '1.5.0',
   default_options : [
     'buildtype=release',
     'warning_level=3',

src/bed.cpp

@@ -0,0 +1,75 @@
+/*!
+ * @file bed.cpp
+ *
+ * @brief BED reader source file
+ */
+
+#include <iostream>
+#include <memory>
+#include <sstream>
+
+#include "bed.hpp"
+
+namespace racon {
+
+bool BedFile::Deserialize(const std::string& line, BedRecord& record) {
+    if (line.empty()) {
+        return false;
+    }
+    char name_buff[1024];
+    int64_t chrom_start = 0, chrom_end = 0;
+    int32_t n = sscanf(line.c_str(), "%s %ld %ld", name_buff, &chrom_start, &chrom_end);
+    if (n < 3 || chrom_end <= chrom_start) {
+        throw std::runtime_error("Invalid BED line: '" + line + "'");
+    }
+    record = BedRecord(name_buff, chrom_start, chrom_end);
+    return true;
+}
+
+void BedFile::Serialize(std::ostream& os, const BedRecord& record) {
+    os << record.chrom() << " " << record.chrom_start() << " " << record.chrom_end();
+}
+
+std::string BedFile::Serialize(const BedRecord& record) {
+    std::ostringstream oss;
+    Serialize(oss, record);
+    return oss.str();
+}
+
+BedReader::BedReader(const std::string& in_fn)
+    : file_{std::unique_ptr<std::ifstream>(new std::ifstream(in_fn))}
+    , in_(*file_.get())
+{
+}
+
+BedReader::BedReader(std::istream& in)
+    : in_(in)
+{
+}
+
+bool BedReader::GetNext(BedRecord& record) {
+    const bool rv1 = !std::getline(in_, line_).fail();
+    if (!rv1)
+        return false;
+    const bool rv2 = BedFile::Deserialize(line_, record);
+    return rv2;
+}
+
+std::vector<BedRecord> BedReader::ReadAll(const std::string& fn)
+{
+    std::ifstream in{fn};
+    return ReadAll(in);
+}
+
+std::vector<BedRecord> BedReader::ReadAll(std::istream& in)
+{
+    std::vector<BedRecord> records;
+    BedReader reader{in};
+    BedRecord record;
+    while (reader.GetNext(record)) {
+        records.emplace_back(std::move(record));
+    }
+    return records;
+}
+
+}

src/bed.hpp

@@ -0,0 +1,93 @@
+/*!
+ * @file bed.hpp
+ *
+ * @brief BED file containers and parser.
+ */
+
+#pragma once
+
+#include <cstdint>
+#include <fstream>
+#include <memory>
+#include <string>
+#include <vector>
+
+namespace racon {
+
+class BedRecord;
+
+class BedFile {
+public:
+    static bool Deserialize(const std::string& line, BedRecord& record);
+    static void Serialize(std::ostream& os, const BedRecord& record);
+    static std::string Serialize(const BedRecord& record);
+};
+
+/*
+ * \brief BedRecord container.
+ * Note: BED records have 0-based coordinates, and the end coordinate is non-inclusive.
+*/
+class BedRecord {
+public:
+    ~BedRecord() = default;
+
+    BedRecord() = default;
+
+    BedRecord(std::string _chrom, int64_t _chrom_start, int64_t _chrom_end)
+        : chrom_(std::move(_chrom))
+        , chrom_start_(_chrom_start)
+        , chrom_end_(_chrom_end) {}
+
+    const std::string& chrom() const {
+        return chrom_;
+    }
+    int64_t chrom_start() const {
+        return chrom_start_;
+    }
+    int64_t chrom_end() const {
+        return chrom_end_;
+    }
+
+    void chrom(const std::string& val) {
+        chrom_ = val;
+    }
+    void chrom_start(int64_t val) {
+        chrom_start_ = val;
+    }
+    void chrom_end(int64_t val) {
+        chrom_end_ = val;
+    }
+
+    bool operator==(const BedRecord& rhs) const
+    {
+        return chrom_ == rhs.chrom_ && chrom_start_ == rhs.chrom_start_
+                    && chrom_end_ == rhs.chrom_end_;
+    }
+
+    std::ostream& operator<<(std::ostream& os) const {
+        BedFile::Serialize(os, *this);
+        return os;
+    }
+
+private:
+    std::string chrom_;
+    int64_t chrom_start_ = 0;
+    int64_t chrom_end_ = 0;
+};
+
+class BedReader {
+public:
+    BedReader(const std::string& in_fn);
+    BedReader(std::istream& in);
+
+    static std::vector<BedRecord> ReadAll(const std::string& fn);
+    static std::vector<BedRecord> ReadAll(std::istream& in);
+    bool GetNext(BedRecord& record);
+
+private:
+    std::unique_ptr<std::ifstream> file_;
+    std::istream& in_;
+    std::string line_;
+};
+
+}

src/main.cpp

@@ -3,6 +3,7 @@
 #include <stdint.h>
 #include <getopt.h>
 
+#include <iostream>
 #include <string>
 #include <vector>
 
@@ -31,6 +32,7 @@ static struct option options[] = {
     {"mismatch", required_argument, 0, 'x'},
     {"gap", required_argument, 0, 'g'},
     {"threads", required_argument, 0, 't'},
+    {"bed", required_argument, 0, 'B'},
     {"version", no_argument, 0, 'v'},
     {"help", no_argument, 0, 'h'},
 #ifdef CUDA_ENABLED
@@ -66,7 +68,9 @@ int main(int argc, char** argv) {
     uint32_t cudaaligner_band_width = 0;
     bool cuda_banded_alignment = false;
 
-    std::string optstring = "ufw:q:e:m:x:g:t:h";
+    std::string bed_file;
+
+    std::string optstring = "ufw:q:e:m:x:g:t:B:h";
 #ifdef CUDA_ENABLED
     optstring += "bc::";
 #endif
@@ -104,6 +108,9 @@ int main(int argc, char** argv) {
             case 't':
                 num_threads = atoi(optarg);
                 break;
+            case 'B':
+                bed_file = std::string(optarg);
+                break;
             case 'v':
                 printf("%s\n", version);
                 exit(0);
@@ -149,8 +156,10 @@ int main(int argc, char** argv) {
         exit(1);
     }
 
+    std::cerr << "BED file: '" << bed_file << "'\n";
+
     auto polisher = racon::createPolisher(input_paths[0], input_paths[1],
-        input_paths[2], type == 0 ? racon::PolisherType::kC :
+        input_paths[2], bed_file, type == 0 ? racon::PolisherType::kC :
         racon::PolisherType::kF, window_length, quality_threshold,
         error_threshold, trim, match, mismatch, gap, num_threads,
         cudapoa_batches, cuda_banded_alignment, cudaaligner_batches,
@@ -209,6 +218,9 @@ void help() {
         "        -g, --gap <int>\n"
         "            default: -4\n"
         "            gap penalty (must be negative)\n"
+        "        -B, --bed <str>\n"
+        "            default: ''\n"
+        "            path to a BED file with regions to polish\n"
         "        -t, --threads <int>\n"
         "            default: 1\n"
         "            number of threads\n"

src/meson.build

@@ -1,9 +1,11 @@
 racon_cpp_sources = files([
+  'bed.cpp',
   'logger.cpp',
   'overlap.cpp',
   'polisher.cpp',
   'sequence.cpp',
-  'window.cpp'
+  'util.cpp',
+  'window.cpp',
 ])
 
 racon_extra_flags = []