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

Author: OndrejSladky

Makefile

@@ -4,13 +4,14 @@ CXX=         g++
 CXXFLAGS=    -g -Wall -Wno-unused-function -std=c++17 -O2 -I/opt/homebrew/include
 LDFLAGS=     -lz -lglpk -L/opt/homebrew/lib
 LARGEFLAGS=  -DLARGE_KMERS
+ELARGEFLAGS=  -DEXTRA_LARGE_KMERS
 SRC=         src
 TESTS=       tests
 GTEST=       $(TESTS)/googletest/googletest
 DATA=        data
 
 
-all: kmercamel kmercamel-large
+all: kmercamel
 
 test: cpptest converttest verify
 
@@ -22,9 +23,10 @@ quick-verify: verify.py kmercamel
 	./verify.py --quick $(DATA)/spneumoniae.fa
 	./verify.py --k 13 --superstring_path $(DATA)/global-k13c.fa $(DATA)/spneumoniae.fa
 
-cpptest: kmercameltest kmercameltest-large
+cpptest: kmercameltest kmercameltest-large kmercameltest-extra-large
 	./kmercameltest
 	./kmercameltest-large
+	./kmercameltest-extra-large
 
 converttest: convert_superstring_unittest.py
 	./convert_superstring_unittest.py
@@ -34,16 +36,15 @@ kmercamel: $(SRC)/main.cpp $(SRC)/$(wildcard *.cpp *.h *.hpp) src/version.h
 	$(CXX) $(CXXFLAGS) $(SRC)/main.cpp -o $@ $(LDFLAGS)
 	cp kmercamel  🐫 || true
 
-kmercamel-large: $(SRC)/main.cpp $(SRC)/$(wildcard *.cpp *.h *.hpp) src/version.h
-	./create-version.sh
-	$(CXX) $(CXXFLAGS) $(SRC)/main.cpp -o $@ $(LDFLAGS) $(LARGEFLAGS)
-
 kmercameltest: $(TESTS)/unittest.cpp gtest-all.o $(SRC)/$(wildcard *.cpp *.h *.hpp) $(TESTS)/$(wildcard *.cpp *.h *.hpp)
 	$(CXX) $(CXXFLAGS) -isystem $(GTEST)/include -I $(GTEST)/include $(TESTS)/unittest.cpp gtest-all.o -pthread -o $@ $(LDFLAGS)
 
 kmercameltest-large: $(TESTS)/unittest.cpp gtest-all.o $(SRC)/$(wildcard *.cpp *.h *.hpp) $(TESTS)/$(wildcard *.cpp *.h *.hpp)
 	$(CXX) $(CXXFLAGS) -isystem $(GTEST)/include -I $(GTEST)/include $(TESTS)/unittest.cpp gtest-all.o -pthread -o $@ $(LDFLAGS) $(LARGEFLAGS)
 
+kmercameltest-extra-large: $(TESTS)/unittest.cpp gtest-all.o $(SRC)/$(wildcard *.cpp *.h *.hpp) $(TESTS)/$(wildcard *.cpp *.h *.hpp)
+	$(CXX) $(CXXFLAGS) -isystem $(GTEST)/include -I $(GTEST)/include $(TESTS)/unittest.cpp gtest-all.o -pthread -o $@ $(LDFLAGS) $(ELARGEFLAGS)
+
 gtest-all.o: $(GTEST)/src/gtest-all.cc $(wildcard *.cpp *.h *.hpp)
 	$(CXX) $(CXXFLAGS) -isystem $(GTEST)/include -I $(GTEST)/include -I $(GTEST) -DGTEST_CREATE_SHARED_LIBRARY=1 -c -pthread $(GTEST)/src/gtest-all.cc -o $@
 
@@ -55,6 +56,7 @@ clean:
 	rm -f kmercamel
 	rm -f 🐫 || true
 	rm -f kmercameltest
+	rm -f kmercameltest-large
 	rm -r -f ./bin
 	rm -f gtest-all.o
 	rm -f src/version.h

README.md

@@ -20,8 +20,7 @@ They come in two different implementations (their results may differ due to the
 Note that at this point only the implementations with hash table are optimized and that the Aho-Corasick automaton
 based versions of the algorithms are only experimental.
 
-The hashing based implementations of the default KmerCamel🐫 (`./kmercamel`) support $k$-mer with $k$ at most 31,
-whereas the larger KmerCamel🐫 (`./kmercamel-large`) supports $k$-mers with $k$ at most 63 (at the cost of slight slowdown).
+The hashing based implementations of KmerCamel🐫 support $k$-mer with $k$ at most 127,
 
 All algorithms can be used to either work in the unidirectional model or in the bidirectional model
 (i.e. treat $k$-mer and its reverse complement as the same; in this case either of them appears in the result).
@@ -62,8 +61,8 @@ on macOS.
 
 The program has the following arguments:
 
-- `-p path_to_fasta` - the path to fasta file. This is a required argument.
-- `-k value_of_k` - the size of one k-mer. This is a required argument.
+- `-p path_to_fasta` - the path to fasta file (can be `gzip`ed). This is a required argument.
+- `-k value_of_k` - the size of one k-mer (up to 127). This is a required argument.
 - `-a algorithm` - the algorithm which should be run. Either `global` or `globalAC` for Global Greedy, `local` or `localAC` for Local Greedy.
 The versions with AC use Aho-Corasick automaton. Default `global`.
 - `-o output_path` - the path to output file. If not specified, output is printed to stdout.
@@ -79,27 +78,18 @@ The output contains the resulting superstring - capital letters indicate that at
 For example:
 
 ```
-./kmercamel -p ./spneumoniae.fa -a local -k 12 -d 7 -c
+./kmercamel -p ./spneumoniae.fa -a local -k 31 -d 5 -c
 ```
 
-runs the Local Greedy in the bidirectional model on the streptococcus fasta file with `k=12` and `d=7`.
+runs the Local Greedy in the bidirectional model on the streptococcus fasta file with `k=31` and `d=5`.
 
 Alternatively, if your operating system supports it, you can run `./🐫` instead of `./kmercamel`.
 
-Currently, KmerCamel🐫 does not support gziped files as an input.
-A possible workaround is to use `gzcat` and process substitution.
-
-```
-./kmercamel -k 13 -p <(gzcat fasta_file.fa.gz)
-```
-
-Note: on some systems you might need to use the name `zcat` instead of `gzcat`.
-
 ### Mask optimization
 
 For mask optimization, run the subcommand `optimize` with the following arguments:
 
-- `p path_to_fasta` - the path to fasta file. This is a required argument.
+- `p path_to_fasta` - the path to fasta file (can be `gzip`ed). 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, `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.
@@ -110,26 +100,11 @@ For mask optimization, run the subcommand `optimize` with the following argument
 For example:
 
 ```
-./kmercamel optimize -p ./global-spneumoniae.fa -k 12 -a runs -c
+./kmercamel optimize -p ./global-spneumoniae.fa -k 31 -a runs -c
 ```
 
 minimizes the number of runs of 1s in the mask of the superstring computed by Global Greedy in the bidirectional model on the streptococcus fasta file with `k=12`.
 
-Note: currently mask optimization cannot read gziped files, nor can the proccess substititution be used.
-
-### Large $k$-mers
-
-The default version of KmerCamel🐫 does not support $k > 31$. For those values, use the large KmerCamel🐫,
-which supports $k < 64$.
-
-For example:
-
-```
-./kmercamel-large -p ./spneumoniae.fa -a global -k 63 -c
-```
-
-Note: for smaller $k$ it is recommended to use default KmerCamel🐫 as it is faster.
-
 ### Turn off memory optimizations for Global
 
 In order to reduce the memory footprint of hash-table based Global Greedy,

src/ac_automaton.h renamed to src/ac/ac_automaton.h

@@ -3,7 +3,7 @@
 #include <list>
 #include <vector>
 
-#include "models.h"
+#include "kmers_ac.h"
 
 constexpr int INVALID_STATE = -1;
 struct ACState {

src/global_ac.h renamed to src/ac/global_ac.h

@@ -10,8 +10,7 @@
 #include <fstream>
 #include <algorithm>
 
-#include "models.h"
-#include "kmers.h"
+#include "kmers_ac.h"
 #include "ac_automaton.h"
 
 

src/ac/kmers_ac.h

@@ -0,0 +1,55 @@
+#pragma once
+#include <string>
+#include <vector>
+#include "../kmers.h"
+
+struct KMer {
+    std::string value;
+    size_t length() const{ return value.size(); }
+};
+
+/// Return the complementary nucleotide for the given one.
+char ComplementaryNucleotide(const char nucleotide) {
+    if (nucleotide == 'A') return 'T';
+    else if (nucleotide == 'T') return 'A';
+    else if (nucleotide == 'G') return 'C';
+    else if (nucleotide == 'C') return 'G';
+    throw std::invalid_argument("cannot find complementary nucleotide for letter " + std::string(1, nucleotide));
+}
+
+/// Compute the reverse complement of the given k-mer.
+KMer ReverseComplement(const KMer &kMer) {
+    KMer ans;
+    ans.value.reserve(kMer.length());
+    for (int i = (int)kMer.length() - 1; i >= 0; --i) {
+        ans.value += ComplementaryNucleotide(kMer.value[i]);
+    }
+    return ans;
+}
+
+/// Convert the given basic nucleotide to int so it can be used for indexing in AC.
+/// If non-existing nucleotide is given, return -1.
+int NucleotideToInt (char c) {
+    switch (c) {
+        case 'A': return 0;
+        case 'C': return 1;
+        case 'G': return 2;
+        case 'T': return 3;
+        case 'a': return 0;
+        case 'c': return 1;
+        case 'g': return 2;
+        case 't': return 3;
+        default: return -1;
+    }
+}
+
+/// Convert the given k-mer to its representation as a number.
+kmer64_t KMerToNumber(const KMer &kMer) {
+    kmer64_t ret = 0;
+    for (char c : kMer.value) {
+        ret <<= 2;
+        ret |= NucleotideToInt(c);
+    }
+    return ret;
+}
+

src/local_ac.h renamed to src/ac/local_ac.h

@@ -7,7 +7,7 @@
 #include <algorithm>
 #include <fstream>
 
-#include "models.h"
+#include "kmers_ac.h"
 #include "ac_automaton.h"
 
 /// Find the index of the first extending k-mer from the incidentKMers which is not forbidden.

src/ac/parser_ac.h

@@ -0,0 +1,79 @@
+#pragma once
+
+// Implicitly initialize kseq.
+#include "../parser.h"
+
+#include "kmers_ac.h"
+
+/// Record of one fasta sequence.
+struct FastaRecord {
+    // The name of the record. Starts with '>'.
+    std::string name;
+    // The genetic sequence.
+    std::string sequence;
+};
+
+/// Read fasta file with given path.
+std::vector<FastaRecord> ReadFasta(std::string &path) {
+    gzFile fp;
+    kseq_t *seq;
+    std::vector<FastaRecord> records;
+
+    fp = gzopen(path.c_str(), "r");
+    if (fp == nullptr) {
+        throw std::invalid_argument("couldn't open file " + path);
+    }
+    seq = kseq_init(fp);
+    while (kseq_read(seq) >= 0) {
+        records.push_back(FastaRecord{
+                seq->name.s,
+                seq->seq.s,
+        });
+    }
+    kseq_destroy(seq);
+    gzclose(fp);
+    return records;
+}
+
+/// Create a list of unique k-mers in no particular order.
+/// This runs in O(k*data.size) expected time.
+void AddKMersFromSequence(std::unordered_set<std::string> &kMers, std::string &data, int k) {
+    // Convert the sequence to uppercase letters.
+    std::transform(data.begin(), data.end(), data.begin(), toupper);
+    size_t possibleKMerEnd = k;
+    for (size_t i = 1; i <= data.size(); ++i) {
+        if (data[i-1] != 'A' && data[i-1] != 'C' && data[i-1] != 'G' && data[i-1] != 'T') {
+            // Skip this and the next k-1 k-mers.
+            possibleKMerEnd = i + k;
+        }
+        if (i >= possibleKMerEnd) {
+            kMers.insert(data.substr(i - k, k));
+        }
+    }
+}
+
+/// Return only the subset of k-mers where no two k-mers are complements of one another.
+/// The k-mers are chosen with no particular property.
+std::unordered_set<std::string> FilterKMersWithComplement(std::unordered_set<std::string> &kMers) {
+    std::unordered_set<std::string> ret;
+    for (auto &&x : kMers) {
+        auto kMer = KMer{x};
+        if (ret.count(ReverseComplement(kMer).value) == 0) ret.insert(x);
+    }
+    return ret;
+}
+
+/// Create a list of unique k-mers in no particular order.
+/// This runs in O(k*data.size) expected time.
+std::vector<KMer> ConstructKMers(std::vector<FastaRecord> &data, int k, bool complements) {
+    std::unordered_set<std::string> uniqueKMers;
+    for (auto &&record : data) {
+        AddKMersFromSequence(uniqueKMers, record.sequence, k);
+    }
+    if (complements) uniqueKMers = FilterKMersWithComplement(uniqueKMers);
+    std::vector<KMer> result;
+    for (auto &kMer : uniqueKMers) {
+        result.push_back(KMer{kMer});
+    }
+    return result;
+}

src/streaming.h renamed to src/ac/streaming.h

@@ -5,7 +5,7 @@
 #include <cassert>
 #include <algorithm>
 
-#include "kmers.h"
+#include "kmers_ac.h"
 
 void Push(std::ostream &of, const std::string &current, int k, int32_t used) {
     if (current.size() == static_cast<size_t>(k) && used != 0) {
@@ -40,7 +40,7 @@ void Streaming(std::string &path, std::ostream &of, int k, bool complements) {
                         if (kMer.size() == size_t(k + 1)) kMer=kMer.substr(1);
                         assert(kMer.size() <= size_t(k));;
                         if (kMer.length() == size_t (k) ) {
-                            int64_t encoded = KMerToNumber(KMer{kMer});
+                            uint64_t encoded = KMerToNumber(KMer{kMer});
                             auto rc = ReverseComplement(encoded, k);
                             auto rc2 = NumberToKMer(rc, k);
                             bool contained = (kMers.count(encoded) > 0) || ((complements) && kMers.count(ReverseComplement(encoded, k)) > 0);