Merge branch 'main' into revert-67-revert-59-lowerbound

Author: PavelVesely

README.md

@@ -101,7 +101,7 @@ For mask optimization, run the subcommand `optimize` with the following argument
 
 - `p path_to_fasta` - the path to fasta file. This is a required argument.
 - `k k_value` - the size of one k-mer. This is a required argument.
-- `a algorithm` - the algorithm for mask optimization. Either `ones` for maximizing the number of 1s, `runs` for minimizing the number of runs of 1s, or `zeros` for maximizing the number of 0s. Default `ones`.
+- `a algorithm` - the algorithm for mask optimization. Either `ones` for maximizing the number of 1s, `runs` for minimizing the number of runs of 1s, `runsapprox` for approximately minimizing the number of runs of 1s, or `zeros` for maximizing the number of 0s. Default `ones`.
 - `o output_path` - the path to output file. If not specified, output is printed to stdout.
 - `c` - treat k-mer and its reverse complement as equal.
 - `h` - print help.

src/main.cpp

@@ -49,7 +49,7 @@ int Help() {
     std::cerr << "Accepted arguments:" << std::endl;
     std::cerr << "  -p path_to_fasta - required; valid path to fasta file" << std::endl;
     std::cerr << "  -k k_value       - required; integer value for k" << std::endl;
-    std::cerr << "  -a algorithm     - the algorithm to be run [ones (default), runs, zeros]" << std::endl;
+    std::cerr << "  -a algorithm     - the algorithm to be run [ones (default), runs, runsapprox, zeros]" << std::endl;
     std::cerr << "  -o output_path   - if not specified, the output is printed to stdout" << std::endl;
     std::cerr << "  -c               - treat k-mer and its reverse complement as equal" << std::endl;
     std::cerr << "  -h               - print help" << std::endl;

src/masks.h

@@ -112,15 +112,15 @@ std::pair<std::vector<int>, std::vector<int>> HeuristicPreSolve(std::vector<std:
 
 
 /// For the given masked superstring output the same superstring with mask with minimal number of runs of ones.
-void OptimizeRuns(std::string path, kh_S64_t *kMers, std::ostream &of, int k, bool complements) {
+void OptimizeRuns(std::string path, kh_S64_t *kMers, std::ostream &of, int k, bool complements, bool approximate) {
     kh_O64_t *intervals = kh_init_O64();
     std::vector<std::list<size_t>> intervalsForKMer;
     auto [size, rows] = ReadIntervals(intervals, kMers, intervalsForKMer, path, k, complements, of, nullptr);
     int mappedSize, newIntervals; size_t totalIntervals;
     auto [mapping, intervalMapping] = HeuristicPreSolve(intervalsForKMer, rows, mappedSize, totalIntervals, newIntervals);
     glp_prob *lp;
     lp = glp_create_prob();
-    if (mappedSize != 0) {
+    if (mappedSize != 0 && !approximate) {
         auto *ia = new int[totalIntervals + 1];
         auto *ja = new int[totalIntervals + 1];
         auto *ar = new double[totalIntervals + 1];
@@ -160,10 +160,10 @@ void OptimizeRuns(std::string path, kh_S64_t *kMers, std::ostream &of, int k, bo
 
     for (size_t i = 0; i < rows; ++i) {
         if (intervalMapping[i] == -1) intervalsSet[i] = true;
+        else if (approximate) intervalsSet[i] = mappedSize != 0;
         else intervalsSet[i] = mappedSize == 0 ? false : (glp_get_col_prim(lp, intervalMapping[i] + 1) > 0.5);
     }
 
-    of << "> superstring" << std::endl;
     ReadIntervals(nullptr, kMers, intervalsForKMer, path, k, complements, of, intervalsSet);
     of << std::endl;
 }
@@ -179,7 +179,9 @@ int Optimize(std::string &algorithm, std::string path, std::ostream &of,  int k,
         } else if (algorithm == "zeros") {
             OptimizeOnes(in, of, kMers, k, complements, true);
         } else if (algorithm == "runs") {
-            OptimizeRuns(path, kMers, of, k, complements);
+            OptimizeRuns(path, kMers, of, k, complements, false);
+        } else if (algorithm == "runsapprox") {
+            OptimizeRuns(path, kMers, of, k, complements, true);
         } else {
             std::cerr << "Algorithm '" + algorithm + "' not recognized." << std::endl;
             in.close();

src/parser.h

@@ -168,7 +168,9 @@ std::pair<size_t, size_t> ReadIntervals(kh_O64_t *intervals, kh_S64_t *kMers, st
                 currentInterval += interval_used;
                 interval_used = false;
             }
-            else if (c == '\n') readingHeader = false;
+            // Reprint the header.
+            if (readingHeader && !reading) of << c;
+            if (c == '\n') readingHeader = false;
             if (readingHeader) continue;
             auto data = NucleotideToInt(c);
             // Disregard white space.
@@ -186,7 +188,7 @@ std::pair<size_t, size_t> ReadIntervals(kh_O64_t *intervals, kh_S64_t *kMers, st
             currentKMer |= data;
             --beforeKMerEnd;
             if (beforeKMerEnd == 0) {
-                bool represented = kh_get_S64(kMers, currentKMer) != kh_end(kMers);
+                bool represented = containsKMer(kMers, currentKMer, k, complements);
                 bool set = false;
                 if (represented) {
                     interval_used = true;

tests/masks_unittest.h

@@ -67,12 +67,13 @@ namespace {
 
     TEST(Mask, OptimizeRuns) {
         std::string path = std::filesystem::current_path();
-        path += "/tests/testdata/runstest.fa";
 
         struct TestCase {
             std::vector<kmer_t> kMers;
             int k;
             bool complements;
+            bool approximate;
+            std::string relativePath;
             std::string wantResult;
         };
         std::vector<TestCase> tests = {
@@ -84,18 +85,42 @@ namespace {
                             KMerToNumber({"TT"}),
                             KMerToNumber({"AT"})
                         },
-                        2, false,
+                        2, false, false,
+                        "/tests/testdata/runstest.fa",
                         "> superstring\nacgcgttACGtATt\n"
                 },
+                {
+                        {
+                                KMerToNumber({"AC"}),
+                                KMerToNumber({"CG"}),
+                                KMerToNumber({"GT"}),
+                                KMerToNumber({"TT"}),
+                                KMerToNumber({"AT"})
+                        },
+                        2, false, true,
+                        "/tests/testdata/runstest.fa",
+                        "> superstring\nACgCGTtACGtATt\n"
+                },
+                {
+                        {
+                                KMerToNumber({"ACT"}),
+                                KMerToNumber({"CTA"}),
+                                KMerToNumber({"GTA"})
+                        },
+                        3, true, false,
+                        "/tests/testdata/runstest2.fa",
+                        "> superstring\nACTAGta\n"
+                },
         };
 
         for (auto &t : tests) {
             std::stringstream of;
+            auto totalPath = path + t.relativePath;
             auto kMersDict = kh_init_S64();
             int ret;
             for (auto &kMer : t.kMers) kh_put_S64(kMersDict, kMer, &ret);
 
-            OptimizeRuns(path, kMersDict, of, t.k, t.complements);
+            OptimizeRuns(totalPath, kMersDict, of, t.k, t.complements, t.approximate);
 
             EXPECT_EQ(t.wantResult, of.str());
         }

tests/testdata/runstest2.fa

@@ -0,0 +1,2 @@
+> superstring
+ACtaGta

verify.py

@@ -84,7 +84,7 @@ def main():
     else:
         k = int(args.k)
         # Do the tests on mask optimization.
-        for a in ["runs", "ones", "zeros"]:
+        for a in ["runs", "ones", "zeros", "runsapprox"]:
             print(f"Testing {a}:")
             for large in [True, False]:
                 if not large and k >= 32: continue