Revert "Revert "KmerCamel can compute lower bounds on superstring lengths.""

Author: OndrejSladky

src/global.h

@@ -57,7 +57,7 @@ void PartialPreSort(std::vector<kmer_t> &vals, int k) {
 /// If complements are provided, treat k-mer and its complement as identical.
 /// If this is the case, k-mers are expected to contain only one k-mer from a complement pair.
 /// Moreover, if so, the resulting Hamiltonian path contains two superstrings which are reverse complements of one another.
-overlapPath OverlapHamiltonianPath (std::vector<kmer_t> &kMers, int k, bool complements) {
+overlapPath OverlapHamiltonianPath (std::vector<kmer_t> &kMers, int k, bool complements, bool lower_bound = false) {
     size_t n = kMers.size();
     size_t kMersCount = n * (1 + complements);
     size_t batchSize = kMersCount / MEMORY_REDUCTION_FACTOR + 1;
@@ -106,10 +106,12 @@ overlapPath OverlapHamiltonianPath (std::vector<kmer_t> &kMers, int k, bool comp
                     if (suffix_key == kh_end(prefixes)) continue;
                     size_t previous, j;
                     previous = j = kh_val(prefixes, suffix_key);
-                    // If the path forms a cycle, or is between k-mer and its reverse complement, or the k-mers complement was already selected skip this path.
                     while (j != size_t(-1) && \
-                           (accessFirstLast(first, last, i, n) % n == j % n \
-                           || accessFirstLast(first, last, i, n) % n == accessFirstLast(last, first, j, n) % n \
+                            // k-mers are complementary
+                           ((i + n) % (2 * n) == j \
+                           // forms a cycle
+                           || (!lower_bound && accessFirstLast(first, last, i, n) == j) \
+                           // k-mer is already used
                            || prefixForbidden[j])) {
                         size_t new_j = next[j - from];
                         // If the k-mer is forbidden, remove it to keep the complexity linear.

src/lower_bound.h

@@ -0,0 +1,14 @@
+#pragma once
+
+#include "global.h"
+#include "kmers.h"
+
+/// Return the length of the cycle cover which lower bounds the superstring length.
+size_t LowerBoundLength(std::vector<kmer_t> kMers, int k, bool complements) {
+    auto cycle_cover = OverlapHamiltonianPath(kMers, k, complements, true);
+    size_t res = 0;
+    for (auto &overlap : cycle_cover.second) {
+        res += size_t(k) - size_t(overlap);
+    }
+    return res / (1 + complements);
+}

src/main.cpp

@@ -11,7 +11,13 @@
 #include "streaming.h"
 #include "output.h"
 #include "khash_utils.h"
+
+#include <iostream>
+#include <string>
+#include "unistd.h"
+#include "version.h"
 #include "masks.h"
+#include "lower_bound.h"
 
 
 #ifdef LARGE_KMERS
@@ -33,6 +39,7 @@ int Help() {
     std::cerr << "  -d d_value       - integer value for d_max; default 5" << std::endl;
     std::cerr << "  -c               - treat k-mer and its reverse complement as equal" << std::endl;
     std::cerr << "  -m               - turn off the memory optimizations for global" << std::endl;
+    std::cerr << "  -l               - compute the cycle cover lower bound instead of masked superstring" << std::endl;
     std::cerr << "  -h               - print help" << std::endl;
     std::cerr << "  -v               - print version" << std::endl;
     std::cerr << "Example usage:       ./kmercamel -p path_to_fasta -k 13 -d 5 -a local" << std::endl;
@@ -71,9 +78,10 @@ int main(int argc, char **argv) {
     bool complements = false;
     bool optimize_memory = true;
     bool d_set = false;
+    bool lower_bound = false;
     int opt;
     try {
-        while ((opt = getopt(argc, argv, "p:k:d:a:o:hcvm"))  != -1) {
+        while ((opt = getopt(argc, argv, "p:k:d:a:o:hcvml"))  != -1) {
             switch(opt) {
                 case  'p':
                     path = optarg;
@@ -103,6 +111,9 @@ int main(int argc, char **argv) {
                 case 'm':
                     optimize_memory = false;
                     break;
+                case 'l':
+                    lower_bound = true;
+                    break;
                 case 'v':
                     Version();
                     return 0;
@@ -140,6 +151,9 @@ int main(int argc, char **argv) {
     } else if (masks && (d_set || !optimize_memory)) {
         std::cerr << "Not supported flags for optimize." << std::endl;
         return Help();
+    } else if (lower_bound && algorithm != "global") {
+        std::cerr << "Lower bound computation supported only for hash table global." << std::endl;
+        return Help();
     }
 
     if (masks) {
@@ -162,14 +176,15 @@ int main(int argc, char **argv) {
             return Help();
         }
         d_max = std::min(k - 1, d_max);
-        WriteName(k, *of);
+        if (!lower_bound) WriteName(k, *of);
         if (algorithm == "global") {
             auto kMerVec = kMersToVec(kMers);
             kh_destroy_S64(kMers);
             // Turn off the memory optimizations if optimize_memory is set to false.
-            if(optimize_memory) PartialPreSort(kMerVec, k);
+            if (optimize_memory) PartialPreSort(kMerVec, k);
             else MEMORY_REDUCTION_FACTOR = 1;
-            Global(kMerVec, *of, k, complements);
+            if (lower_bound) std::cout << LowerBoundLength(kMerVec, k, complements);
+            else Global(kMerVec, *of, k, complements);
         }
         else Local(kMers, *of, k, d_max, complements);
     } else {

tests/global_unittest.h

@@ -78,30 +78,48 @@ namespace {
             overlapPath wantResult;
             int k;
             bool complements;
+            bool lower_bound;
         };
         std::vector<TestCase> tests = {
                 {
                         {KMerToNumber({"AT"})},
                         {{(size_t)-1, (size_t)-1}, {(byte)-1, (byte)-1}},
                         2,
                         true,
+                        false,
                 },
                 {
                         {KMerToNumber({"ACG"}), KMerToNumber({"TAC"}), KMerToNumber({"GGC"})},
                         {{2, 0, (size_t)-1}, {1, 2, (byte)-1}},
                         3,
                         false,
+                        false,
                 },
                 {
                         {KMerToNumber({"ACAA"}), KMerToNumber({"ATTT"}), KMerToNumber({"AACA"})},
                         {{4, 3, 0, 5, (size_t)-1, (size_t)-1}, {2, 2, 3, 3, (byte)-1, (byte)-1}},
                         4,
                         true,
+                        false,
+                },
+                {
+                        {KMerToNumber({"ACG"}), KMerToNumber({"CGT"}), KMerToNumber({"TAA"})},
+                        {{1, 2, 0}, {2, 1, 1}},
+                        3,
+                        false,
+                        true,
+                },
+                {
+                        {KMerToNumber({"ACC"}), KMerToNumber({"CGG"})},
+                        {{3, 2, 1, 0}, {2, 2, 0, 0}},
+                        3,
+                        true,
+                        true,
                 },
         };
 
         for (auto t : tests) {
-            overlapPath got = OverlapHamiltonianPath( t.kMers, t.k, t.complements);
+            overlapPath got = OverlapHamiltonianPath( t.kMers, t.k, t.complements, t.lower_bound);
             EXPECT_EQ(t.wantResult.first, got.first);
             EXPECT_EQ(t.wantResult.second, got.second);
         }

tests/lower_bound_unittest.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include "../src/lower_bound.h"
+
+#include "gtest/gtest.h"
+
+namespace {
+    TEST(LowerBound, LowerBoundLength) {
+        struct TestCase {
+            std::vector<kmer_t> kMers;
+            int k;
+            bool complements;
+            size_t wantResult;
+        };
+        std::vector<TestCase> tests = {
+                {
+                        {KMerToNumber({"ACG"}), KMerToNumber({"CGT"}), KMerToNumber({"TAA"})},
+                        3, false, 5
+                },
+                {
+                        {KMerToNumber({"AC"}), KMerToNumber({"GG"})},
+                        2, true, 3
+                },
+                {
+                        {KMerToNumber({"AAACCC"}), KMerToNumber({"CCCAAA"}), KMerToNumber({"TGGGGT"})},
+                        6, false, 11
+                },
+        };
+
+        for (auto &t : tests) {
+            auto gotResult = LowerBoundLength(t.kMers, t.k, t.complements);
+
+            ASSERT_EQ(t.wantResult, gotResult);
+        }
+    }
+}

tests/unittest.cpp

@@ -6,6 +6,7 @@
 #include "local_ac_unittest.h"
 #include "streaming_unittest.h"
 #include "ac_automaton_unittest.h"
+#include "lower_bound_unittest.h"
 #include "masks_unittest.h"
 
 #include "gtest/gtest.h"