[FEATURE] Allow writing and reading a sketch file (#254)

* [FEATURE] Read and write a sketch file in chopper layout.

* [MISC] automatic linting

* Review

* [MISC] automatic linting

---------

Co-authored-by: seqan-actions[bot] <[email protected]>
Co-authored-by: Enrico Seiler <[email protected]>

Author: 3 people

Diff for: include/chopper/sketch/read_hll_files_into.hpp

@@ -0,0 +1,45 @@
+// ---------------------------------------------------------------------------------------------------
+// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin
+// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik
+// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
+// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md
+// ---------------------------------------------------------------------------------------------------
+
+#pragma once
+
+#include <vector>
+
+#include <cereal/types/vector.hpp>
+
+#include <chopper/configuration.hpp>
+
+#include <hibf/sketch/hyperloglog.hpp>
+#include <hibf/sketch/minhashes.hpp>
+
+namespace chopper::sketch
+{
+
+struct sketch_file
+{
+    chopper::configuration chopper_config{};
+    std::vector<std::vector<std::string>> filenames{};
+    std::vector<seqan::hibf::sketch::hyperloglog> hll_sketches{};
+    std::vector<seqan::hibf::sketch::minhashes> minHash_sketches{};
+
+private:
+    friend class cereal::access;
+
+    template <typename archive_t>
+    void serialize(archive_t & archive)
+    {
+        uint32_t version{1};
+        archive(CEREAL_NVP(version));
+
+        archive(CEREAL_NVP(chopper_config));
+        archive(CEREAL_NVP(filenames));
+        archive(CEREAL_NVP(hll_sketches));
+        archive(CEREAL_NVP(minHash_sketches));
+    }
+};
+
+} // namespace chopper::sketch

Diff for: include/chopper/sketch/sketch_file.hpp

@@ -21,12 +21,31 @@
 #include <chopper/sketch/check_filenames.hpp>
 #include <chopper/sketch/output.hpp>
 #include <chopper/sketch/read_data_file.hpp>
+#include <chopper/sketch/sketch_file.hpp>
 
 #include <hibf/sketch/compute_sketches.hpp>
 
 namespace chopper
 {
 
+void validate_configuration(sharg::parser & parser,
+                            chopper::configuration & config,
+                            chopper::configuration const & sketch_config)
+{
+    if (parser.is_option_set("sketch-bits"))
+        throw sharg::parser_error{"You cannot set --sketch-bits when using a sketch file as input."};
+
+    if (parser.is_option_set("kmer") && config.k != sketch_config.k)
+    {
+        std::cerr << sharg::detail::to_string(
+            "[WARNING] Given k-mer size (",
+            config.k,
+            ") differs from k-mer size in the sketch file (",
+            config.k,
+            "). The results may be suboptimal. If this was a conscious decision, you can ignore this warning.\n");
+    }
+}
+
 int chopper_layout(chopper::configuration & config, sharg::parser & parser)
 {
     parser.parse();
@@ -36,41 +55,70 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser)
     else if (config.k > config.window_size)
         throw sharg::parser_error{"The k-mer size cannot be bigger than the window size."};
 
+    auto has_sketch_file_extension = [](std::filesystem::path const & path)
+    {
+        return path.string().ends_with(".sketch") || path.string().ends_with(".sketches");
+    };
+
     config.disable_sketch_output = !parser.is_option_set("output-sketches-to");
+    if (!config.disable_sketch_output && !has_sketch_file_extension(config.sketch_directory))
+        throw sharg::parser_error{"The sketch output file must have the extension \".sketch\" or \".sketches\"."};
+
+    bool const input_is_a_sketch_file = has_sketch_file_extension(config.data_file);
 
     int exit_code{};
 
     std::vector<std::vector<std::string>> filenames{};
+    std::vector<seqan::hibf::sketch::hyperloglog> sketches{};
 
-    chopper::sketch::read_data_file(config, filenames);
+    if (input_is_a_sketch_file)
+    {
+        chopper::sketch::sketch_file sin{};
+
+        { // Deserialization is guaranteed to be complete when going out of scope.
+            std::ifstream is{config.data_file};
+            cereal::BinaryInputArchive iarchive{is};
+            iarchive(sin);
+        }
 
-    std::vector<seqan::hibf::sketch::hyperloglog> sketches;
+        filenames = std::move(sin.filenames); // No need to call check_filenames because the files are not read.
+        sketches = std::move(sin.hll_sketches);
+        validate_configuration(parser, config, sin.chopper_config);
+    }
+    else
+    {
+        chopper::sketch::read_data_file(config, filenames);
 
-    if (filenames.empty())
-        throw sharg::parser_error{
-            sharg::detail::to_string("The file ", config.data_file.string(), " appears to be empty.")};
+        if (filenames.empty())
+            throw sharg::parser_error{
+                sharg::detail::to_string("The file ", config.data_file.string(), " appears to be empty.")};
 
-    chopper::sketch::check_filenames(filenames, config);
+        // Files need to exist because they will be read for sketching.
+        chopper::sketch::check_filenames(filenames, config);
+    }
 
     config.hibf_config.input_fn =
         chopper::input_functor{filenames, config.precomputed_files, config.k, config.window_size};
     config.hibf_config.number_of_user_bins = filenames.size();
     config.hibf_config.validate_and_set_defaults();
 
-    config.compute_sketches_timer.start();
-    seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches);
-    config.compute_sketches_timer.stop();
+    if (!input_is_a_sketch_file)
+    {
+        config.compute_sketches_timer.start();
+        seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches);
+        config.compute_sketches_timer.stop();
+    }
 
     exit_code |= chopper::layout::execute(config, filenames, sketches);
 
     if (!config.disable_sketch_output)
     {
-        if (!std::filesystem::exists(config.sketch_directory))
-            std::filesystem::create_directory(config.sketch_directory);
-
-        assert(filenames.size() == sketches.size());
-        for (size_t i = 0; i < filenames.size(); ++i)
-            chopper::sketch::write_sketch_file(filenames[i][0], sketches[i], config);
+        chopper::sketch::sketch_file sout{.chopper_config = config,
+                                          .filenames = std::move(filenames),
+                                          .hll_sketches = std::move(sketches)};
+        std::ofstream os{config.sketch_directory, std::ios::binary};
+        cereal::BinaryOutputArchive oarchive{os};
+        oarchive(sout);
     }
 
     if (!config.output_timings.empty())

Diff for: src/chopper_layout.cpp

@@ -42,9 +42,10 @@ void set_up_parser(sharg::parser & parser, configuration & config)
             .short_id = '\0',
             .long_id = "input",
             .description =
-                "The input must be a file containing paths to sequence data you wish to estimate; one filepath "
-                "per line. If your file contains auxiliary information (e.g. species IDs), your file must be tab-"
-                "separated.",
+                "The input can either be (1) a sketch (binary) file (extensions: \".sketch\" or \".sketches\") that "
+                "was produced by chopper or (2) a file containing paths to sequence data you wish to estimate; one "
+                "filepath per line. If your file contains auxiliary information (e.g. species IDs), your file must be "
+                "tab-separated.",
             .required = true});
     parser.add_list_item("", "Example file:");
     parser.add_list_item("", "```");
@@ -248,8 +249,8 @@ void set_up_parser(sharg::parser & parser, configuration & config)
         sharg::config{
             .long_id = "output-sketches-to",
             .description =
-                "If you supply a directory path with this option, the hyperloglog sketches of your input will be "
-                "stored in the respective path; one .hll file per input file.",
+                "If supplied, a sketch file will be written to the specified directory. The sketch file can be used "
+                "as input for the next run of chopper. The file extension must be \".sketch\" or \".sketches\".",
             .default_message = "None",
             .advanced = true});
 

Diff for: src/set_up_parser.cpp

@@ -1,7 +1,7 @@
 cmake_minimum_required (VERSION 3.18)
 
 if (NOT TARGET chopper_sketch)
-    add_library (chopper_sketch STATIC check_filenames.cpp output.cpp read_data_file.cpp read_hll_files_into.cpp)
+    add_library (chopper_sketch STATIC check_filenames.cpp output.cpp read_data_file.cpp)
 
     target_link_libraries (chopper_sketch PUBLIC chopper_shared)
 endif ()

Diff for: src/sketch/CMakeLists.txt

@@ -21,7 +21,6 @@
 
 #include <chopper/configuration.hpp>
 #include <chopper/layout/execute.hpp>
-#include <chopper/sketch/read_hll_files_into.hpp>
 
 #include <hibf/sketch/compute_sketches.hpp>
 #include <hibf/sketch/hyperloglog.hpp>

Diff for: src/sketch/read_hll_files_into.cpp

@@ -7,12 +7,6 @@
 
 cmake_minimum_required (VERSION 3.18)
 
-add_api_test (read_hll_files_into_test.cpp)
-target_use_datasources (read_hll_files_into_test FILES small.fa)
-target_use_datasources (read_hll_files_into_test FILES small2.fa)
-target_use_datasources (read_hll_files_into_test FILES small.hll)
-target_use_datasources (read_hll_files_into_test FILES small2.hll)
-
 add_api_test (check_filenames_test.cpp)
 target_use_datasources (check_filenames_test FILES seq1.fa)
 target_use_datasources (check_filenames_test FILES seq2.fa)