Giraffe integration for haplotype sampling

src/recombinator.cpp

@@ -1033,15 +1033,46 @@ double get_or_estimate_coverage(
             count_to_frequency[iter->second]++;
         }
     }
+
+    // Use mode as the initial estimate for coverage.
     auto statistics = summary_statistics(count_to_frequency);
-    double coverage = std::max(statistics.median, statistics.mode);
+    double coverage = statistics.mode;
+    bool reliable = true;
     if (verbosity >= Haplotypes::verbosity_detailed) {
         std::cerr << "Coverage: median " << statistics.median
             << ", mean " << statistics.mean
             << ", stdev " << statistics.stdev
-            << ", mode " << statistics.mode
-            << "; using " << coverage << std::endl;
+            << ", mode " << statistics.mode;
+    }
+
+    // If mode < median, try to find a secondary peak at ~2x mode and use
+    // it if it is good enough.
+    if (statistics.mode < statistics.median) {
+        size_t low = 1.7 * statistics.mode, high = 2.3 * statistics.mode;
+        size_t peak = count_to_frequency[coverage];
+        size_t best = low, secondary = count_to_frequency[low];
+        for (size_t i = low + 1; i <= high; i++) {
+            if (count_to_frequency[i] > secondary) {
+                best = i; secondary = count_to_frequency[i];
+            }
+        }
+        if (verbosity >= Haplotypes::verbosity_detailed) {
+            std::cerr << "; secondary peak at " << best;
+        }
+        if (best >= size_t(statistics.median) && secondary >= peak / 2) {
+            coverage = best;
+        } else {
+            reliable = false;
+        }
     }
+
+    if (verbosity >= Haplotypes::verbosity_detailed) {
+        std::cerr << "; using " << coverage << std::endl;
+    }
+    if (!reliable) {
+        std::cerr << "warning: Kmer coverage estimate is unreliable" << std::endl;
+    }
+
     if (verbosity >= Haplotypes::verbosity_basic) {
         double seconds = gbwt::readTimer() - start;
         std::cerr << "Estimated kmer coverage in " << seconds << " seconds" << std::endl;

src/subcommand/haplotypes_main.cpp

@@ -389,7 +389,7 @@ HaplotypesConfig::HaplotypesConfig(int argc, char** argv, size_t max_threads) {
             {
                 std::string arg = optarg;
                 size_t offset = arg.find(':');
-                if (offset == 0 || offset + 1 >= arg.length()) {
+                if (offset == 0 || offset == std::string::npos || offset + 1 >= arg.length()) {
                     std::cerr << "error: [vg haplotypes] cannot parse chain:subchain from " << arg << std::endl;
                     std::exit(EXIT_FAILURE);
                 }
@@ -519,6 +519,7 @@ void preprocess_graph(const gbwtgraph::GBZ& gbz, Haplotypes& haplotypes, Haploty
         validate_haplotypes(haplotypes, gbz.graph, r_index, minimizer_index, expected_chains, config.verbosity);
     }
 }
+
 //----------------------------------------------------------------------------
 
 size_t threads_to_jobs(size_t threads) {

src/utility.cpp

@@ -648,6 +648,15 @@ pair<string, string> split_ext(const string& filename) {
     return parts;
 }
 
+string file_base_name(const string& filename) {
+    size_t slash = filename.rfind('/');
+    if (slash == string::npos) {
+        return split_ext(filename).first;
+    } else {
+        return split_ext(filename.substr(slash + 1)).first;
+    }
+}
+
 bool file_exists(const string& filename) {
     // TODO: use C++17 features to actually poll existence.
     // For now we see if we can open it.

src/utility.hpp

@@ -808,6 +808,9 @@ void get_input_file(const string& file_name, function<void(istream&)> callback);
 /// Split off the extension from a filename and return both parts. 
 pair<string, string> split_ext(const string& filename);
 
+/// Get the base name of a filename (without the directory and the extension).
+string file_base_name(const string& filename);
+
 /// Determine if a file exists.
 /// Only works for files readable by the current user.
 bool file_exists(const string& filename);