spoa-upgrade: added min_coverage argument

README.md

@@ -95,6 +95,9 @@ Usage of `racon` is as following:
             maximum allowed error rate used for filtering overlaps
         --no-trimming
             disables consensus trimming at window ends
+        --min-coverage <int>
+            default: -1
+            minimal consensus coverage
         -m, --match <int>
             default: 3
             score for matching bases

src/main.cpp

@@ -22,6 +22,7 @@ static struct option options[] = {
     {"quality-threshold", required_argument, 0, 'q'},
     {"error-threshold", required_argument, 0, 'e'},
     {"no-trimming", no_argument, 0, 'T'},
+    {"min-coverage", required_argument, nullptr, 'M'},
     {"match", required_argument, 0, 'm'},
     {"mismatch", required_argument, 0, 'x'},
     {"gap", required_argument, 0, 'g'},
@@ -48,6 +49,7 @@ int main(int argc, char** argv) {
     double error_threshold = 0.3;
     bool trim = true;
 
+    int32_t min_coverage = -1;
     int8_t match = 3;
     int8_t mismatch = -5;
     int8_t gap = -4;
@@ -87,6 +89,9 @@ int main(int argc, char** argv) {
             case 'T':
                 trim = false;
                 break;
+            case 'M':
+                min_coverage = atoi(optarg);
+                break;
             case 'm':
                 match = atoi(optarg);
                 break;
@@ -147,7 +152,7 @@ int main(int argc, char** argv) {
     auto polisher = racon::createPolisher(input_paths[0], input_paths[1],
         input_paths[2], type == 0 ? racon::PolisherType::kC :
         racon::PolisherType::kF, window_length, quality_threshold,
-        error_threshold, trim, match, mismatch, gap, num_threads,
+        error_threshold, trim, min_coverage, match, mismatch, gap, num_threads,
         cudapoa_batches, cuda_banded_alignment, cudaaligner_batches,
         cudaaligner_band_width);
 
@@ -195,6 +200,9 @@ void help() {
         "            maximum allowed error rate used for filtering overlaps\n"
         "        --no-trimming\n"
         "            disables consensus trimming at window ends\n"
+        "        --min-coverage <int>\n"
+        "            default: -1\n"
+        "            minimal consensus coverage\n"
         "        -m, --match <int>\n"
         "            default: 3\n"
         "            score for matching bases\n"

src/polisher.cpp

@@ -58,8 +58,9 @@ std::unique_ptr<Polisher> createPolisher(const std::string& sequences_path,
     const std::string& overlaps_path, const std::string& target_path,
     PolisherType type, uint32_t window_length, double quality_threshold,
     double error_threshold, bool trim, int8_t match, int8_t mismatch, int8_t gap,
-    uint32_t num_threads, uint32_t cudapoa_batches, bool cuda_banded_alignment,
-    uint32_t cudaaligner_batches, uint32_t cudaaligner_band_width) {
+    uint32_t num_threads, int32_t min_coverage, uint32_t cudapoa_batches,
+    bool cuda_banded_alignment, uint32_t cudaaligner_batches,
+    uint32_t cudaaligner_band_width) {
 
     if (type != PolisherType::kC && type != PolisherType::kF) {
         fprintf(stderr, "[racon::createPolisher] error: invalid polisher type!\n");
@@ -158,7 +159,7 @@ std::unique_ptr<Polisher> createPolisher(const std::string& sequences_path,
         return std::unique_ptr<Polisher>(new Polisher(std::move(sparser),
                     std::move(oparser), std::move(tparser), type, window_length,
                     quality_threshold, error_threshold, trim, match, mismatch, gap,
-                    num_threads));
+                    num_threads, min_coverage));
     }
 }
 
@@ -167,13 +168,13 @@ Polisher::Polisher(std::unique_ptr<bioparser::Parser<Sequence>> sparser,
     std::unique_ptr<bioparser::Parser<Sequence>> tparser,
     PolisherType type, uint32_t window_length, double quality_threshold,
     double error_threshold, bool trim, int8_t match, int8_t mismatch, int8_t gap,
-    uint32_t num_threads)
+    uint32_t num_threads, int32_t min_coverage)
         : sparser_(std::move(sparser)), oparser_(std::move(oparser)),
         tparser_(std::move(tparser)), type_(type), quality_threshold_(
         quality_threshold), error_threshold_(error_threshold), trim_(trim),
-        alignment_engines_(), sequences_(), dummy_quality_(window_length, '!'),
-        window_length_(window_length), windows_(),
-        thread_pool_(std::make_shared<thread_pool::ThreadPool>(num_threads)),
+        min_coverage_(min_coverage), alignment_engines_(), sequences_(),
+        dummy_quality_(window_length, '!'), window_length_(window_length),
+        windows_(), thread_pool_(std::make_shared<thread_pool::ThreadPool>(num_threads)),
         logger_(new Logger()) {
 
     for (uint32_t i = 0; i < num_threads; ++i) {
@@ -498,7 +499,7 @@ void Polisher::polish(std::vector<std::unique_ptr<Sequence>>& dst,
             [&](uint64_t j) -> bool {
                 auto it = thread_pool_->thread_map().find(std::this_thread::get_id());  // NOLINT
                 return windows_[j]->generate_consensus(
-                    alignment_engines_[it->second], trim_);
+                    alignment_engines_[it->second], trim_, min_coverage_);
             }, i));
     }
 

src/polisher.hpp

@@ -44,7 +44,7 @@ std::unique_ptr<Polisher> createPolisher(const std::string& sequences_path,
     const std::string& overlaps_path, const std::string& target_path,
     PolisherType type, uint32_t window_length, double quality_threshold,
     double error_threshold, bool trim, int8_t match, int8_t mismatch, int8_t gap,
-    uint32_t num_threads, uint32_t cuda_batches = 0,
+    uint32_t num_threads, int32_t min_coverage = -1, uint32_t cuda_batches = 0,
     bool cuda_banded_alignment = false, uint32_t cudaaligner_batches = 0,
     uint32_t cudaaligner_band_width = 0);
 
@@ -61,16 +61,17 @@ class Polisher {
         const std::string& overlaps_path, const std::string& target_path,
         PolisherType type, uint32_t window_length, double quality_threshold,
         double error_threshold, bool trim, int8_t match, int8_t mismatch, int8_t gap,
-        uint32_t num_threads, uint32_t cuda_batches, bool cuda_banded_alignment,
-        uint32_t cudaaligner_batches, uint32_t cudaaligner_band_width);
+        uint32_t num_threads, int32_t min_coverage, uint32_t cuda_batches,
+        bool cuda_banded_alignment, uint32_t cudaaligner_batches,
+        uint32_t cudaaligner_band_width);
 
 protected:
     Polisher(std::unique_ptr<bioparser::Parser<Sequence>> sparser,
         std::unique_ptr<bioparser::Parser<Overlap>> oparser,
         std::unique_ptr<bioparser::Parser<Sequence>> tparser,
         PolisherType type, uint32_t window_length, double quality_threshold,
         double error_threshold, bool trim, int8_t match, int8_t mismatch, int8_t gap,
-        uint32_t num_threads);
+        uint32_t num_threads, int32_t min_coverage);
     Polisher(const Polisher&) = delete;
     const Polisher& operator=(const Polisher&) = delete;
     virtual void find_overlap_breaking_points(std::vector<std::unique_ptr<Overlap>>& overlaps);
@@ -83,6 +84,7 @@ class Polisher {
     double quality_threshold_;
     double error_threshold_;
     bool trim_;
+    int32_t min_coverage_;
     std::vector<std::shared_ptr<spoa::AlignmentEngine>> alignment_engines_;
 
     std::vector<std::unique_ptr<Sequence>> sequences_;

src/window.cpp

@@ -63,7 +63,7 @@ void Window::add_layer(const char* sequence, uint32_t sequence_length,
 }
 
 bool Window::generate_consensus(std::shared_ptr<spoa::AlignmentEngine> alignment_engine,
-    bool trim) {
+    bool trim, int32_t min_coverage) {
 
     if (sequences_.size() < 3) {
         consensus_ = std::string(sequences_.front().first, sequences_.front().second);
@@ -120,7 +120,7 @@ bool Window::generate_consensus(std::shared_ptr<spoa::AlignmentEngine> alignment
     }
 
     std::vector<uint32_t> coverages;
-    consensus_ = graph.GenerateConsensus(&coverages);
+    consensus_ = graph.GenerateConsensus(&coverages, min_coverage);
 
     if (type_ == WindowType::kTGS && trim) {
         uint32_t average_coverage = (sequences_.size() - 1) / 2;

src/window.hpp

@@ -45,7 +45,7 @@ class Window {
     }
 
     bool generate_consensus(std::shared_ptr<spoa::AlignmentEngine> alignment_engine,
-        bool trim);
+        bool trim, int32_t min_coverage);
 
     void add_layer(const char* sequence, uint32_t sequence_length,
         const char* quality, uint32_t quality_length, uint32_t begin,