diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5c59b3b..5a60f8f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.11)
 project(raven 
-          VERSION 1.8.1
+          VERSION 1.8.2
           LANGUAGES CXX
           DESCRIPTION "Raven is a de novo genome assembler for long uncorrected reads.")
 
diff --git a/Raven.deps.cmake b/Raven.deps.cmake
index e32b1b0..10c656b 100644
--- a/Raven.deps.cmake
+++ b/Raven.deps.cmake
@@ -1,6 +1,20 @@
 include(FetchContent)
 include(GNUInstallDirs)
 
+FetchContent_Declare(
+  ram
+  GIT_REPOSITORY https://github.com/tbrekalo/ram-experimental
+  GIT_TAG f7042c21fee01feb87e1953e498c402eb1f8fcc6)
+
+FetchContent_GetProperties(ram)
+if (NOT ram_POPULATED)
+  FetchContent_Populate(ram)
+  add_subdirectory(
+    ${ram_SOURCE_DIR} 
+    ${ram_BINARY_DIR}
+    EXCLUDE_FROM_ALL)
+endif ()
+
 find_package(bioparser 3.0.13 QUIET)
 if (NOT bioparser_FOUND)
     FetchContent_Declare(
@@ -53,6 +67,40 @@ if (NOT racon_FOUND)
     endif ()
 endif ()
 
+find_package(cxxopts 3.0.0 QUIET)
+if (NOT cxxopts_FOUND)
+  FetchContent_Declare(
+    cxxopts
+    GIT_REPOSITORY https://github.com/jarro2783/cxxopts.git
+    GIT_TAG v3.0.0)
+
+  FetchContent_GetProperties(cxxopts)
+  if (NOT cxxopts_POPULATED)
+    FetchContent_Populate(cxxopts)
+    add_subdirectory(
+      ${cxxopts_SOURCE_DIR}
+      ${cxxopts_BINARY_DIR}
+      EXCLUDE_FROM_ALL)
+  endif ()
+endif ()
+
+find_package(tsl-robin-map 1.0.1 QUIET)
+if (NOT tsl-robin-map_FOUND)
+  FetchContent_Declare(
+    tsl-robin-map
+    GIT_REPOSITORY https://github.com/Tessil/robin-map.git
+    GIT_TAG v1.0.1)
+
+  FetchContent_GetProperties(tsl-robin-map)
+  if (NOT tsl-robin-map_POPULATED)
+    FetchContent_Populate(tsl-robin-map)
+    add_subdirectory(
+      ${tsl-robin-map_SOURCE_DIR}
+      ${tsl-robin-map_BINARY_DIR}
+      EXCLUDE_FROM_ALL)
+  endif ()
+endif ()
+
 if (RAVEN_BUILD_PYTHON)
     find_package(pybind11 QUIET)
     if (NOT pybind11_FOUND)
diff --git a/RavenExe/Exe.cmake b/RavenExe/Exe.cmake
index fa588b6..178077d 100644
--- a/RavenExe/Exe.cmake
+++ b/RavenExe/Exe.cmake
@@ -3,7 +3,7 @@ if (NOT TARGET RavenExe)
     include(${CMAKE_CURRENT_LIST_DIR}/Exe.srcs.cmake)
 
     add_executable(raven_exe ${SOURCES})
-    target_link_libraries(raven_exe PRIVATE raven)
+    target_link_libraries(raven_exe PRIVATE cxxopts::cxxopts raven)
 
     if (racon_enable_cuda)
         target_compile_definitions(raven_exe PRIVATE CUDA_ENABLED)
@@ -12,4 +12,7 @@ if (NOT TARGET RavenExe)
     target_compile_definitions(raven_exe PRIVATE VERSION="${PROJECT_VERSION}")
     set_property(TARGET raven_exe PROPERTY OUTPUT_NAME raven)
 
+    configure_file(
+      "${CMAKE_CURRENT_LIST_DIR}/src/raven_cfg.h.in" raven_cfg.h)
+    target_include_directories(raven_exe PRIVATE ${PROJECT_BINARY_DIR})
 endif()
diff --git a/RavenExe/src/main.cc b/RavenExe/src/main.cc
index 7472013..d8649b3 100644
--- a/RavenExe/src/main.cc
+++ b/RavenExe/src/main.cc
@@ -1,213 +1,119 @@
 // Copyright (c) 2020 Robert Vaser
 
-#include <getopt.h>
-
-#include <iostream>
 #include <cstdlib>
+#include <filesystem>
+#include <iostream>
+#include <string_view>
 
 #include "bioparser/fasta_parser.hpp"
 #include "bioparser/fastq_parser.hpp"
 #include "biosoup/timer.hpp"
+#include "cxxopts.hpp"
 #include "raven/raven.h"
-std::atomic<std::uint32_t> biosoup::NucleicAcid::num_objects{0};
-
-namespace {
-
-static struct option options[] = {
-    {"kmer-len", required_argument, nullptr, 'k'},
-    {"window-len", required_argument, nullptr, 'w'},
-    {"frequency", required_argument, nullptr, 'f'},
-    {"identity", required_argument, nullptr, 'i'},
-    {"kMaxNumOverlaps", required_argument, nullptr, 'o'},
-    {"polishing-rounds", required_argument, nullptr, 'p'},
-    {"match", required_argument, nullptr, 'm'},
-    {"mismatch", required_argument, nullptr, 'n'},
-    {"gap", required_argument, nullptr, 'g'},
-#ifdef CUDA_ENABLED
-    {"cuda-poa-batches", optional_argument, nullptr, 'c'},
-    {"cuda-banded-alignment", no_argument, nullptr, 'b'},
-    {"cuda-alignment-batches", required_argument, nullptr, 'a'},
-#endif
-    {"graphical-fragment-assembly", required_argument, nullptr, 'F'},
-    {"resume", no_argument, nullptr, 'r'},
-    {"disable-checkpoints", no_argument, nullptr, 'd'},
-    {"threads", required_argument, nullptr, 't'},
-    {"version", no_argument, nullptr, 'v'},
-    {"help", no_argument, nullptr, 'h'},
-    {nullptr, 0, nullptr, 0}};
-
-void Help() {
-  std::cout
-      << "usage: raven [options ...] <sequences> [<sequences> ...]\n"
-         "\n"
-         "  # default output is to stdout in FASTA format\n"
-         "  <sequences>\n"
-         "    input file in FASTA/FASTQ format (can be compressed with gzip)\n"
-         "\n"
-         "  options:\n"
-         "    -k, --kmer-len <int>\n"
-         "      default: 15\n"
-         "      length of minimizers used to find overlaps\n"
-         "    -w, --window-len <int>\n"
-         "      default: 5\n"
-         "      length of sliding window from which minimizers are sampled\n"
-         "    -f, --frequency <double>\n"
-         "      default: 0.001\n"
-         "      threshold for ignoring most frequent minimizers\n"
-         "    -i, --identity <double>\n"
-         "      default: 0\n"
-         "      threshold for overlap between two reads in order to construct an edge between them\n"
-         "    -o, --kMaxNumOverlaps <long unsigned int>\n"
-         "      default: 32\n"
-         "      maximum number of overlaps that will be taken during FindOverlapsAndCreatePiles stage\n"
-         "    -p, --polishing-rounds <int>\n"
-         "      default: 2\n"
-         "      number of times racon is invoked\n"
-         "    -m, --match <int>\n"
-         "      default: 3\n"
-         "      score for matching bases\n"
-         "    -n, --mismatch <int>\n"
-         "      default: -5\n"
-         "      score for mismatching bases\n"
-         "    -g, --gap <int>\n"
-         "      default: -4\n"
-         "      gap penalty (must be negative)\n"
-#ifdef CUDA_ENABLED
-         "    -c, --cuda-poa-batches <int>\n"
-         "      default: 0\n"
-         "      number of batches for CUDA accelerated polishing\n"
-         "    -b, --cuda-banded-alignment\n"
-         "      use banding approximation for polishing on GPU\n"
-         "      (only applicable when -c is used)\n"
-         "    -a, --cuda-alignment-batches <int>\n"
-         "      default: 0\n"
-         "      number of batches for CUDA accelerated alignment\n"
-#endif
-         "    --graphical-fragment-assembly <string>\n"
-         "      prints the assembly graph in GFA format\n"
-         "    --resume\n"
-         "      resume previous run from last checkpoint\n"
-         "    --disable-checkpoints\n"
-         "      disable checkpoint file creation\n"
-         "    -t, --threads <int>\n"
-         "      default: 1\n"
-         "      number of threads\n"
-         "    --version\n"
-         "      prints the version number\n"
-         "    -h, --help\n"
-         "      prints the usage\n";
-}
+#include "raven_cfg.h"
 
-}  // namespace
+std::atomic<std::uint32_t> biosoup::NucleicAcid::num_objects{0};
 
 int main(int argc, char** argv) {
-  std::uint8_t kmer_len = 15;
-  std::uint8_t window_len = 5;
-  double freq = 0.001;
-  double identity = 0;
-  std::size_t kMaxNumOverlaps = 32;
-
-  std::uint32_t num_polishing_rounds = 2;
-  std::int8_t m = 3;
-  std::int8_t n = -5;
-  std::int8_t g = -4;
+  auto options = cxxopts::Options(
+      "raven",
+      "Raven is a de novo genome assembler for long uncorrected reads");
+
+  /* clang-format off */
+  options.add_options("maping arguments")
+    ("k,kmer-len", "length of minimizers used to find overlaps",
+      cxxopts::value<std::uint8_t>()->default_value("15"))
+    ("w,window-len", "length of sliding window from which minimizers are sampled",
+      cxxopts::value<std::uint8_t>()->default_value("5"))
+    ("f,frequency", "threshold for ignoring most frequent minimizers",
+      cxxopts::value<double>()->default_value("0.0002"));
+  options.add_options("layout arguments")
+    ("i,identity", 
+     "threshold for overlap between two reads in order to construct an edge between them", 
+      cxxopts::value<double>()->default_value("0"))
+    ("o,kMaxNumOverlaps", 
+      "maximum number of overlaps that will be taken during FindOverlapsAndCreatePiles stage", 
+      cxxopts::value<std::size_t>()->default_value("32"));
+  options.add_options("polishing arguments")
+    ("p,polishing-rounds", "number of times racon is invoked",
+      cxxopts::value<std::uint32_t>()->default_value("2"))
+    ("m,match", "score for matching bases",
+      cxxopts::value<std::int8_t>()->default_value("3"))
+    ("n,mismatch", "score for mismatching bases",
+      cxxopts::value<std::int8_t>()->default_value("-5"))
+    ("g,gap", "gap penalty (must be negative)",
+      cxxopts::value<std::int8_t>()->default_value("-4"));
+
+#ifndef CUDA_ENABLED
+  options.add_options("cuda arguments")
+    ("c,cuda-poa-batches", "number of batches for CUDA accelerated polishing",
+      cxxopts::value<std::uint32_t>()->default_value("0"))
+    ("b,cuda-banded-alignment", 
+     "use banding approximation for polishing on GPU; (only applicable when -c is used)",
+     cxxopts::value<bool>()->default_value("false"))
+    ("a,cuda-alignment-batches", "number of batches for CUDA accelerated alignment",
+      cxxopts::value<std::uint32_t>()->default_value("0"));
+#endif /* CUDA_ENABLED */
+
+  options.add_options("utility arguments")
+    ("graphical-fragment-assembly", "prints the assembly graph in GFA format",
+      cxxopts::value<std::string>()->default_value("./raven.gfa"))
+    ("disable-checkpoints", "disable checkpoint file creation",
+      cxxopts::value<bool>()->default_value("false"))
+    ("resume", "resume previous run from last checkpoint",
+      cxxopts::value<bool>()->default_value("false"))
+    ("version", "prints the version number",
+      cxxopts::value<bool>()->default_value("false"))
+    ("t,threads", "number of threads", 
+      cxxopts::value<std::uint32_t>()->default_value("1"))
+    ("h,help", "display help message");
+
+  options.add_options("sequence input")
+    ("paths", "input fasta/fastq reads", 
+      cxxopts::value<std::vector<std::string>>());
+  /* clang-format on */
+
+  options.parse_positional("paths");
+  auto result = options.parse(argc, argv);
+
+  if (result.count("help")) {
+    std::cerr << options.help() << std::endl;
+    return EXIT_SUCCESS;
+  }
 
-  std::string gfa_path = "";
-  bool resume = false;
-  bool checkpoints = true;
+  if (result.count("version")) {
+    std::cerr << PROJECT_VER << std::endl;
+    return EXIT_SUCCESS;
+  }
 
-  std::uint32_t num_threads = 1;
+  const auto kmer_len = result["kmer-len"].as<std::uint8_t>();
+  const auto window_len = result["window-len"].as<std::uint8_t>();
+  const auto freq = result["frequency"].as<double>();
+  const auto identity = result["identity"].as<double>();
+  const auto kMaxNumOverlaps = result["kMaxNumOverlaps"].as<std::size_t>();
 
-  std::uint32_t cuda_poa_batches = 0;
-  std::uint32_t cuda_alignment_batches = 0;
-  bool cuda_banded_alignment = false;
+  const auto num_polishing_rounds =
+      result["polishing-rounds"].as<std::uint32_t>();
+  const auto m = result["match"].as<std::int8_t>();
+  const auto n = result["mismatch"].as<std::int8_t>();
+  const auto g = result["gap"].as<std::int8_t>();
 
-  std::string optstr = "k:w:f:p:m:n:g:t:h";
-#ifdef CUDA_ENABLED
-  optstr += "c:ba:";
-#endif
-  int arg;
-  while ((arg = getopt_long(argc, argv, optstr.c_str(), options, nullptr)) !=
-         -1) {  // NOLINT
-    switch (arg) {
-      case 'k':
-        kmer_len = std::atoi(optarg);
-        break;
-      case 'w':
-        window_len = std::atoi(optarg);
-        break;
-      case 'f':
-        freq = std::atof(optarg);
-        break;
-      case 'i':
-        identity = std::atof(optarg);
-        break;
-      case 'o':
-        kMaxNumOverlaps = std::atof(optarg);
-        break;
-      case 'p':
-        num_polishing_rounds = std::atoi(optarg);
-        break;
-      case 'm':
-        m = std::atoi(optarg);
-        break;
-      case 'n':
-        n = std::atoi(optarg);
-        break;
-      case 'g':
-        g = std::atoi(optarg);
-        break;
-#ifdef CUDA_ENABLED
-      case 'c':
-        cuda_poa_batches = 1;
-        // next text entry is not an option, assuming it's the arg for option c
-        if (optarg == NULL && argv[optind] != NULL && argv[optind][0] != '-') {
-          cuda_poa_batches = std::atoi(argv[optind++]);
-        }
-        // optional argument provided in the ususal way
-        if (optarg != NULL) {
-          cuda_poa_batches = std::atoi(optarg);
-        }
-        break;
-      case 'b':
-        cuda_banded_alignment = true;
-        break;
-      case 'a':
-        cuda_alignment_batches = std::atoi(optarg);
-        break;
-#endif
-      case 'F':
-        gfa_path = optarg;
-        break;
-      case 'r':
-        resume = true;
-        break;
-      case 'd':
-        checkpoints = false;
-        break;
-      case 't':
-        num_threads = std::atoi(optarg);
-        break;
-      case 'v':
-        std::cout << VERSION << std::endl;
-        return EXIT_SUCCESS;
-      case 'h':
-        Help();
-        return EXIT_SUCCESS;
-      default:
-        return EXIT_FAILURE;
-    }
-  }
+  const auto gfa_path = result["graphical-fragment-assembly"].as<std::string>();
+  const auto resume = result["resume"].as<bool>();
+  const auto checkpoints = !result["disable-checkpoints"].as<bool>();
 
-  if (argc == 1) {
-    Help();
-    return EXIT_SUCCESS;
-  }
+  const auto num_threads = result["threads"].as<std::uint32_t>();
 
-  if (optind >= argc) {
-    std::cerr << "[raven::] error: missing input file(s)!" << std::endl;
-    return EXIT_FAILURE;
-  }
+#ifdef CUDA_ENABLED
+  const auto cuda_poa_batches = 0U;
+  const auto cuda_alignment_batches = 0U;
+  const auto cuda_banded_alignment = false;
+#else
+  const auto cuda_poa_batches = result["cuda-poa-batches"].as<std::uint32_t>();
+  const auto cuda_alignment_batches =
+      result["cuda-alignment-batches"].as<std::uint32_t>();
+  const auto cuda_banded_alignment = result["cuda-banded-alignment"].as<bool>();
+#endif /* CUDA_ENABLED */
 
   biosoup::Timer timer{};
   timer.Start();
@@ -232,35 +138,25 @@ int main(int argc, char** argv) {
 
   std::vector<std::unique_ptr<biosoup::NucleicAcid>> sequences;
   if (graph.stage < -3 || num_polishing_rounds > std::max(0, graph.stage)) {
-    for (int i = optind; i < argc; ++i) {
-      auto sparser = raven::CreateParser(argv[i]);
-      if (sparser == nullptr) {
-        return EXIT_FAILURE;
-      }
-
-      std::vector<std::unique_ptr<biosoup::NucleicAcid>> chunk;
-      try {
-        chunk = sparser->Parse(-1);
-      } catch (const std::invalid_argument& exception) {
-        std::cerr << exception.what() << " (" << argv[i] << ")" << std::endl;
-        return EXIT_FAILURE;
-      }
-
-      if (chunk.empty()) {
-        std::cerr << "[raven::] warning: file " << argv[i] << " is empty"
-                  << std::endl;
-        continue;
-      }
-
-      if (sequences.empty()) {
-        sequences.swap(chunk);
-      } else {
-        sequences.insert(sequences.end(),
-                         std::make_move_iterator(chunk.begin()),
-                         std::make_move_iterator(chunk.end()));
+    auto paths_strs = result["paths"].as<std::vector<std::string>>();
+    auto paths = std::vector<std::filesystem::path>();
+
+    for (const auto& it : paths_strs) {
+      auto it_path = std::filesystem::path(std::move(it));
+      if (std::filesystem::is_regular_file(it_path)) {
+        paths.emplace_back(std::move(it_path));
+      } else if (std::filesystem::is_directory(it_path)) {
+        for (auto dir_entry : std::filesystem::recursive_directory_iterator(
+                 std::move(it_path))) {
+          if (std::filesystem::is_regular_file(dir_entry)) {
+            paths.emplace_back(std::move(dir_entry));
+          }
+        }
       }
     }
 
+    sequences = raven::LoadSequences(thread_pool, paths);
+
     if (sequences.empty()) {
       std::cerr << "[raven::] error: empty sequences set" << std::endl;
       return EXIT_FAILURE;
@@ -275,8 +171,11 @@ int main(int argc, char** argv) {
 
   raven::ConstructGraph(
       graph, sequences, thread_pool, checkpoints,
-      raven::OverlapPhaseCfg{
-          .kmer_len = kmer_len, .window_len = window_len, .freq = freq, .identity = identity, .kMaxNumOverlaps = kMaxNumOverlaps});
+      raven::OverlapPhaseCfg{.kmer_len = kmer_len,
+                             .window_len = window_len,
+                             .freq = freq,
+                             .identity = identity,
+                             .kMaxNumOverlaps = kMaxNumOverlaps});
 
   raven::Assemble(thread_pool, graph, checkpoints);
 
diff --git a/RavenExe/src/raven_cfg.h.in b/RavenExe/src/raven_cfg.h.in
new file mode 100644
index 0000000..7c2bfea
--- /dev/null
+++ b/RavenExe/src/raven_cfg.h.in
@@ -0,0 +1,10 @@
+#ifndef RAVEN_CFG_H_IN_
+#define RAVEN_CFG_H_IN_
+
+#define PROJECT_NAME "@PROJECT_NAME@"
+#define PROJECT_VER  "@PROJECT_VERSION@"
+#define PROJECT_VER_MAJOR "@PROJECT_VERSION_MAJOR@"
+#define PROJECT_VER_MINOR "@PROJECT_VERSION_MINOR@"
+#define PTOJECT_VER_PATCH "@PROJECT_VERSION_PATCH@"
+
+#endif // RAVEN_CFG_H_IN_
diff --git a/RavenLib/RavenLib.cmake b/RavenLib/RavenLib.cmake
index c1a1f10..75b29e9 100644
--- a/RavenLib/RavenLib.cmake
+++ b/RavenLib/RavenLib.cmake
@@ -37,5 +37,6 @@ if (NOT TARGET RavenLib)
       bioparser::bioparser
       cereal::cereal 
       ram::ram
-      racon::racon)
+      racon::racon
+      robin_map)
 endif ()
diff --git a/RavenLib/include/raven/graph/common.h b/RavenLib/include/raven/graph/common.h
index e658448..c1b36b3 100644
--- a/RavenLib/include/raven/graph/common.h
+++ b/RavenLib/include/raven/graph/common.h
@@ -12,7 +12,7 @@ namespace raven {
 // - ignore nodes less than epsilon away from a node with outdegree > 1
 RAVEN_EXPORT std::uint32_t CreateUnitigs(Graph& graph, std::uint32_t epsilon = 0);
 
-RAVEN_EXPORT void RemoveEdges(Graph& graph, const std::unordered_set<std::uint32_t>& indices,
+RAVEN_EXPORT void RemoveEdges(Graph& graph, const tsl::robin_set<std::uint32_t>& indices,
                  bool remove_nodes = false);
 
 RAVEN_EXPORT std::vector<std::unique_ptr<biosoup::NucleicAcid>> GetUnitigs(
diff --git a/RavenLib/include/raven/graph/graph.h b/RavenLib/include/raven/graph/graph.h
index 81ab86c..f6b47d5 100644
--- a/RavenLib/include/raven/graph/graph.h
+++ b/RavenLib/include/raven/graph/graph.h
@@ -5,20 +5,10 @@
 #include <unordered_set>
 
 #include "biosoup/nucleic_acid.hpp"
+#include "tsl/robin_set.h"
 #include "raven/export.h"
 #include "raven/pile.h"
 
-namespace biosoup {
-
-template <class Archive>
-void serialize(Archive& archive, NucleicAcid& sequence) {
-  archive(sequence.id, sequence.name, sequence.deflated_data,
-          sequence.block_quality, sequence.inflated_len,
-          sequence.is_reverse_complement);
-}
-
-}  // namespace biosoup
-
 namespace raven {
 
 namespace detail {
@@ -118,19 +108,13 @@ struct Node {
     return outdegree() > 0 && indegree() == 0 && count < 6;
   }
 
-  template <class Archive>
-  void serialize(Archive& archive) {
-    archive(id, sequence, count, is_unitig, is_circular, is_polished,
-            transitive);
-  }
-
   std::uint32_t id;
   biosoup::NucleicAcid sequence;
   std::uint32_t count;
   bool is_unitig;
   bool is_circular;
   bool is_polished;
-  std::unordered_set<std::uint32_t> transitive;
+  tsl::robin_set<std::uint32_t> transitive;
   std::vector<Edge*> inedges;
   std::vector<Edge*> outedges;
   Node* pair;
@@ -150,11 +134,6 @@ struct Edge {
 
   bool is_rc() const { return id & 1; }
 
-  template <class Archive>
-  void serialize(Archive& archive) {  // NOLINT
-    archive(id, length, weight);
-  }
-
   std::uint32_t id;
   std::uint32_t length;
   double weight;
diff --git a/RavenLib/include/raven/io.h b/RavenLib/include/raven/io.h
index 6e27010..fcd764c 100644
--- a/RavenLib/include/raven/io.h
+++ b/RavenLib/include/raven/io.h
@@ -1,18 +1,24 @@
 #ifndef RAVEN_IO_H_
 #define RAVEN_IO_H_
 
-#include <string_view>
+#include <filesystem>
+#include <vector>
 
-#include "bioparser/fasta_parser.hpp"
-#include "bioparser/fastq_parser.hpp"
-#include "raven/export.h"
 #include "biosoup/nucleic_acid.hpp"
+#include "raven/export.h"
+#include "thread_pool/thread_pool.hpp"
 
 namespace raven {
 
-RAVEN_EXPORT std::unique_ptr<bioparser::Parser<biosoup::NucleicAcid>> CreateParser(
-    const std::string& path);
+RAVEN_EXPORT
+std::vector<std::unique_ptr<biosoup::NucleicAcid>> LoadSequences(
+    const std::filesystem::path& path);
+
+RAVEN_EXPORT
+std::vector<std::unique_ptr<biosoup::NucleicAcid>> LoadSequences(
+    std::shared_ptr<thread_pool::ThreadPool> thread_pool,
+    const std::vector<std::filesystem::path>& paths);
 
-}
+}  // namespace raven
 
 #endif  // RAVEN_IO_H_
diff --git a/RavenLib/src/assemble.cc b/RavenLib/src/assemble.cc
index c6ba657..99c5686 100644
--- a/RavenLib/src/assemble.cc
+++ b/RavenLib/src/assemble.cc
@@ -33,7 +33,7 @@ static std::uint32_t RemoveTransitiveEdges(Graph& graph) {
   };
 
   std::vector<Edge*> candidate(graph.nodes.size(), nullptr);
-  std::unordered_set<std::uint32_t> marked_edges;
+  tsl::robin_set<std::uint32_t> marked_edges;
 
   for (const auto& it : graph.nodes) {
     if (it == nullptr) {
@@ -104,7 +104,7 @@ static std::uint32_t RemoveTips(Graph& graph) {
       continue;
     }
 
-    std::unordered_set<std::uint32_t> marked_edges;
+    tsl::robin_set<std::uint32_t> marked_edges;
     for (auto jt : end->outedges) {
       if (jt->head->indegree() > 1) {
         marked_edges.emplace(jt->id);
@@ -127,10 +127,10 @@ static std::uint32_t RemoveTips(Graph& graph) {
   return num_tips;
 }
 
-static std::unordered_set<std::uint32_t> FindRemovableEdges(
+static tsl::robin_set<std::uint32_t> FindRemovableEdges(
     const std::vector<Node*>& path) {
   if (path.empty()) {
-    return std::unordered_set<std::uint32_t>{};
+    return tsl::robin_set<std::uint32_t>{};
   }
 
   auto find_edge = [](Node* tail, Node* head) -> Edge* {
@@ -159,7 +159,7 @@ static std::unordered_set<std::uint32_t> FindRemovableEdges(
     }
   }
 
-  std::unordered_set<std::uint32_t> dst;
+  tsl::robin_set<std::uint32_t> dst;
   if (pref == -1 && suff == -1) {  // remove whole path
     for (std::uint32_t i = 0; i < path.size() - 1; ++i) {
       auto it = find_edge(path[i], path[i + 1]);
@@ -239,35 +239,35 @@ static std::uint32_t RemoveBubbles(Graph& graph) {
   };
   auto bubble_pop =
       [&](const std::vector<Node*>& lhs,
-          const std::vector<Node*>& rhs) -> std::unordered_set<std::uint32_t> {
+          const std::vector<Node*>& rhs) -> tsl::robin_set<std::uint32_t> {
     if (lhs.empty() || rhs.empty()) {
-      return std::unordered_set<std::uint32_t>{};
+      return tsl::robin_set<std::uint32_t>{};
     }
 
     // check BFS result
-    std::unordered_set<Node*> bubble;
+    tsl::robin_set<Node*> bubble;
     bubble.insert(lhs.begin(), lhs.end());
     bubble.insert(rhs.begin(), rhs.end());
     if (lhs.size() + rhs.size() - 2 != bubble.size()) {
-      return std::unordered_set<std::uint32_t>{};
+      return tsl::robin_set<std::uint32_t>{};
     }
     for (const auto& it : lhs) {
       if (bubble.count(it->pair) != 0) {
-        return std::unordered_set<std::uint32_t>{};
+        return tsl::robin_set<std::uint32_t>{};
       }
     }
 
     if (!path_is_simple(lhs) || !path_is_simple(rhs)) {  // complex path(s)
       // check poppability
       if (FindRemovableEdges(lhs).empty() && FindRemovableEdges(rhs).empty()) {
-        return std::unordered_set<std::uint32_t>{};
+        return tsl::robin_set<std::uint32_t>{};
       }
 
       // check similarity
       auto l = path_sequence(lhs);
       auto r = path_sequence(rhs);
       if (std::min(l.size(), r.size()) < std::max(l.size(), r.size()) * 0.8) {
-        return std::unordered_set<std::uint32_t>{};
+        return tsl::robin_set<std::uint32_t>{};
       }
 
       EdlibAlignResult result = edlibAlign(l.c_str(), l.size(), r.c_str(),
@@ -279,7 +279,7 @@ static std::uint32_t RemoveBubbles(Graph& graph) {
         edlibFreeAlignResult(result);
       }
       if (score < 0.8) {
-        return std::unordered_set<std::uint32_t>{};
+        return tsl::robin_set<std::uint32_t>{};
       }
     }
 
@@ -336,7 +336,7 @@ static std::uint32_t RemoveBubbles(Graph& graph) {
       }
     }
 
-    std::unordered_set<std::uint32_t> marked_edges;
+    tsl::robin_set<std::uint32_t> marked_edges;
     if (end) {
       auto lhs = path_extract(begin, end);
       auto rhs = path_extract(begin, other_end);
@@ -366,7 +366,7 @@ static void CreateForceDirectedLayout(
     os << "{" << std::endl;
   }
 
-  std::vector<std::unordered_set<std::uint32_t>> components;
+  std::vector<tsl::robin_set<std::uint32_t>> components;
   std::vector<char> is_visited(graph.nodes.size(), 0);
   for (std::uint32_t i = 0; i < graph.nodes.size(); ++i) {
     if (graph.nodes[i] == nullptr || is_visited[i]) {
@@ -399,8 +399,8 @@ static void CreateForceDirectedLayout(
   std::vector<char>().swap(is_visited);
 
   std::sort(components.begin(), components.end(),
-            [](const std::unordered_set<std::uint32_t>& lhs,
-               const std::unordered_set<std::uint32_t>& rhs) {
+            [](const tsl::robin_set<std::uint32_t>& lhs,
+               const tsl::robin_set<std::uint32_t>& rhs) {
               return lhs.size() > rhs.size();
             });
 
@@ -526,7 +526,7 @@ static void CreateForceDirectedLayout(
 
     // update transitive edges
     for (const auto& n : component) {
-      std::unordered_set<std::uint32_t> valid;
+      tsl::robin_set<std::uint32_t> valid;
       for (const auto& m : graph.nodes[n]->transitive) {
         if (component.find(m) != component.end()) {
           valid.emplace(m);
@@ -707,7 +707,7 @@ static std::uint32_t RemoveLongEdges(
   for (std::uint32_t i = 0; i < num_rounds; ++i) {
     CreateForceDirectedLayout(threadPool, graph);
 
-    std::unordered_set<std::uint32_t> marked_edges;
+    tsl::robin_set<std::uint32_t> marked_edges;
     for (const auto& it : graph.nodes) {
       if (it == nullptr || it->outdegree() < 2) {
         continue;
diff --git a/RavenLib/src/binary.cc b/RavenLib/src/binary.cc
index 3e7f64d..f7622f5 100644
--- a/RavenLib/src/binary.cc
+++ b/RavenLib/src/binary.cc
@@ -6,18 +6,39 @@
 #include "cereal/archives/binary.hpp"
 #include "cereal/types/memory.hpp"
 #include "cereal/types/string.hpp"
-#include "cereal/types/unordered_set.hpp"
 #include "cereal/types/vector.hpp"
 
+namespace cereal {
+
+template <class Archive>
+void serialize(Archive& archive, biosoup::NucleicAcid& sequence) {
+  archive(sequence.id, sequence.name, sequence.deflated_data,
+          sequence.block_quality, sequence.inflated_len,
+          sequence.is_reverse_complement);
+}
+
 template <class Archive>
 void serialize(Archive& archive, raven::Edge& edge) {
   archive(edge.id, edge.length, edge.weight);
 }
 
 template <class Archive>
-void serialize(Archive& archive, raven::Node& node) {
+void save(Archive& archive, const raven::Node& node) {
+  auto trans_buff = std::vector<std::uint32_t>(node.transitive.cbegin(),
+                                               node.transitive.cend());
+
   archive(node.id, node.sequence, node.count, node.is_unitig, node.is_circular,
-          node.is_polished, node.transitive);
+          node.is_polished, trans_buff);
+}
+
+template <class Archive>
+void load(Archive& archive, raven::Node& node) {
+  auto trans_buff = std::vector<std::uint32_t>();
+  archive(node.id, node.sequence, node.count, node.is_unitig, node.is_circular,
+          node.is_polished, trans_buff);
+
+  node.transitive =
+      tsl::robin_set<std::uint32_t>(trans_buff.cbegin(), trans_buff.cend());
 }
 
 template <class Archive>
@@ -68,6 +89,8 @@ void load(Archive& archive, raven::Graph& graph) {
   graph.edge_factory = raven::EdgeFactory(graph.edges.size());
 }
 
+}  // namespace cereal
+
 namespace raven {
 
 void StoreGraphToFile(const Graph& graph) {
diff --git a/RavenLib/src/common.cc b/RavenLib/src/common.cc
index 7721b5e..56da9a3 100644
--- a/RavenLib/src/common.cc
+++ b/RavenLib/src/common.cc
@@ -2,13 +2,13 @@
 
 namespace raven {
 
-void RemoveEdges(Graph& graph, const std::unordered_set<std::uint32_t>& indices,
+void RemoveEdges(Graph& graph, const tsl::robin_set<std::uint32_t>& indices,
                  bool remove_nodes) {
   auto erase_remove = [](std::vector<Edge*>& edges, Edge* marked) -> void {
     edges.erase(std::remove(edges.begin(), edges.end(), marked), edges.end());
   };
 
-  std::unordered_set<std::uint32_t> node_indices;
+  tsl::robin_set<std::uint32_t> node_indices;
   for (auto i : indices) {
     if (remove_nodes) {
       node_indices.emplace(graph.edges[i]->tail->id);
@@ -30,7 +30,7 @@ void RemoveEdges(Graph& graph, const std::unordered_set<std::uint32_t>& indices,
 }
 
 std::uint32_t CreateUnitigs(Graph& graph, std::uint32_t epsilon) {
-  std::unordered_set<std::uint32_t> marked_edges;
+  tsl::robin_set<std::uint32_t> marked_edges;
   std::vector<std::unique_ptr<Node>> unitigs;
   std::vector<std::unique_ptr<Edge>> unitig_edges;
   std::vector<std::uint32_t> node_updates(graph.nodes.size(), 0);
@@ -171,7 +171,7 @@ std::uint32_t CreateUnitigs(Graph& graph, std::uint32_t epsilon) {
 
   for (const auto& it : graph.nodes) {  // update transitive edges
     if (it) {
-      std::unordered_set<std::uint32_t> valid;
+      tsl::robin_set<std::uint32_t> valid;
       for (auto jt : it->transitive) {
         valid.emplace(node_updates[jt] == 0 ? jt : node_updates[jt]);
       }
diff --git a/RavenLib/src/construct.cc b/RavenLib/src/construct.cc
index 27cfea2..6b9b8b3 100644
--- a/RavenLib/src/construct.cc
+++ b/RavenLib/src/construct.cc
@@ -1,13 +1,14 @@
 #include "raven/graph/construct.h"
 
+#include <cmath>
 #include <deque>
 
 #include "biosoup/overlap.hpp"
 #include "biosoup/timer.hpp"
+#include "edlib.h"
 #include "ram/minimizer_engine.hpp"
 #include "raven/graph/overlap_utils.h"
 #include "raven/graph/serialization/binary.h"
-#include "edlib.h"
 
 namespace raven {
 
@@ -18,7 +19,6 @@ void FindOverlapsAndCreatePiles(
     double freq, std::vector<std::unique_ptr<Pile>>& piles,
     std::vector<std::vector<biosoup::Overlap>>& overlaps,
     std::size_t kMaxNumOverlaps) {
-  
   piles.reserve(sequences.size());
   for (const auto& it : sequences) {
     piles.emplace_back(new Pile(it->id, it->inflated_len));
@@ -84,26 +84,88 @@ void FindOverlapsAndCreatePiles(
         }
 
         void_futures.emplace_back(thread_pool->Submit(
-            [&](std::uint32_t i) -> void {
-              piles[i]->AddLayers(overlaps[i].begin() + num_overlaps[i],
-                                  overlaps[i].end());
+            [&](std::uint32_t read_id) -> void {
+              piles[read_id]->AddLayers(
+                  overlaps[read_id].begin() + num_overlaps[read_id],
+                  overlaps[read_id].end());
 
-              num_overlaps[i] = std::min(overlaps[i].size(), kMaxNumOverlaps);
+              num_overlaps[read_id] =
+                  std::min(overlaps[read_id].size(), kMaxNumOverlaps);
 
-              if (overlaps[i].size() < kMaxNumOverlaps) {
+              if (overlaps[read_id].size() < kMaxNumOverlaps) {
                 return;
               }
 
-              std::sort(overlaps[i].begin(), overlaps[i].end(),
-                        [&](const biosoup::Overlap& lhs,
-                            const biosoup::Overlap& rhs) -> bool {
-                          return GetOverlapLength(lhs) > GetOverlapLength(rhs);
-                        });
+              auto const kSegOvlspN = static_cast<std::size_t>(
+                  std::round(static_cast<double>(kMaxNumOverlaps) / 3.0));
+
+              auto const kBounds = std::array<std::uint32_t, 4U>{
+                  0U,
+                  static_cast<std::uint32_t>(sequences[read_id]->inflated_len *
+                                             0.333),
+                  static_cast<std::uint32_t>(sequences[read_id]->inflated_len *
+                                             0.666),
+                  sequences[read_id]->inflated_len};
+
+              auto windows = std::array<
+                  std::vector<std::pair<std::uint32_t, std::uint32_t>>, 3U>();
+              for (auto& window : windows) {
+                window.reserve(kSegOvlspN + 1U);
+              }
+
+              auto valid_ids = std::vector<std::uint32_t>();
+              valid_ids.reserve(kMaxNumOverlaps + 3U);
+
+              for (auto ovlp_id = 0U; ovlp_id < overlaps[read_id].size(); ++ovlp_id) {
+                for (auto i = 0U; i < windows.size(); ++i) {
+                  if (overlaps[read_id][ovlp_id].lhs_begin < kBounds[i + 1U] &&
+                      overlaps[read_id][ovlp_id].lhs_end >= kBounds[i]) {
+                    windows[i].emplace_back(
+                        GetOverlapLength(overlaps[read_id][ovlp_id]), ovlp_id);
+                    for (auto j = windows[i].size() - 1U; j > 0; --j) {
+                      if (windows[i][j - 1U].first < windows[i][j].first) {
+                        std::swap(windows[i][j - 1U], windows[i][j]);
+                      } else {
+                        break;
+                      }
+
+                      if (windows[i].size() > kSegOvlspN) {
+                        windows[i].pop_back();
+                      }
+                    }
+                  }
+                }
+              }
+
+              for (auto const& window : windows) {
+                for (auto const& [score, id] : window) {
+                  valid_ids.push_back(id);
+                }
+              }
+
+              std::sort(valid_ids.begin(), valid_ids.end());
+              valid_ids.erase(std::unique(valid_ids.begin(), valid_ids.end()),
+                              valid_ids.end());
+
+              auto updated_ovlps = std::vector<biosoup::Overlap>();
+              updated_ovlps.reserve(valid_ids.size());
 
-              std::vector<biosoup::Overlap> tmp;
-              tmp.insert(tmp.end(), overlaps[i].begin(),
-                         overlaps[i].begin() + kMaxNumOverlaps);  // NOLINT
-              tmp.swap(overlaps[i]);
+              for (auto const id : valid_ids) {
+                updated_ovlps.push_back(overlaps[read_id][id]);
+              }
+
+              std::swap(updated_ovlps, overlaps[read_id]);
+              // std::sort(overlaps[i].begin(), overlaps[i].end(),
+              //           [&](const biosoup::Overlap& lhs,
+              //               const biosoup::Overlap& rhs) -> bool {
+              //             return GetOverlapLength(lhs) >
+              //             GetOverlapLength(rhs);
+              //           });
+
+              // std::vector<biosoup::Overlap> tmp;
+              // tmp.insert(tmp.end(), overlaps[i].begin(),
+              //            overlaps[i].begin() + kMaxNumOverlaps);  // NOLINT
+              // tmp.swap(overlaps[i]);
             },
             it->id()));
       }
@@ -164,7 +226,7 @@ void ResolveContainedReads(
     std::vector<std::future<void>> futures;
     for (std::uint32_t i = 0; i < overlaps.size(); ++i) {
       futures.emplace_back(thread_pool->Submit(
-          [&] (std::uint32_t i) -> void {
+          [&](std::uint32_t i) -> void {
             std::uint32_t k = 0;
             for (std::uint32_t j = 0; j < overlaps[i].size(); ++j) {
               if (!OverlapUpdate(overlaps[i][j], piles)) {
@@ -174,26 +236,24 @@ void ResolveContainedReads(
               const auto& it = overlaps[i][j];
 
               auto lhs = sequences[it.lhs_id]->InflateData(
-                  it.lhs_begin,
-                  it.lhs_end - it.lhs_begin);
+                  it.lhs_begin, it.lhs_end - it.lhs_begin);
 
               auto rhs = sequences[it.rhs_id]->InflateData(
-                  it.rhs_begin,
-                  it.rhs_end - it.rhs_begin);
+                  it.rhs_begin, it.rhs_end - it.rhs_begin);
               if (!it.strand) {
                 biosoup::NucleicAcid rhs_{"", rhs};
                 rhs_.ReverseAndComplement();
                 rhs = rhs_.InflateData();
               }
 
-              auto result = edlibAlign(
-                  lhs.c_str(), lhs.size(),
-                  rhs.c_str(), rhs.size(),
-                  edlibDefaultAlignConfig());
+              auto result = edlibAlign(lhs.c_str(), lhs.size(), rhs.c_str(),
+                                       rhs.size(), edlibDefaultAlignConfig());
 
-              auto score = result.status == EDLIB_STATUS_OK ?
-                  1. - static_cast<double>(result.editDistance) / std::max(lhs.size(), rhs.size()) :  // NOLINT
-                  0.;
+              auto score = result.status == EDLIB_STATUS_OK
+                               ? 1. - static_cast<double>(result.editDistance) /
+                                          std::max(lhs.size(), rhs.size())
+                               :  // NOLINT
+                               0.;
 
               edlibFreeAlignResult(result);
 
@@ -210,9 +270,8 @@ void ResolveContainedReads(
       it.wait();
     }
 
-    std::cerr << "[raven::Graph::Construct] filtered overlaps "
-              << std::fixed << timer.Stop() << "s"
-              << std::endl;
+    std::cerr << "[raven::Graph::Construct] filtered overlaps " << std::fixed
+              << timer.Stop() << "s" << std::endl;
   }
 
   timer.Start();
@@ -314,8 +373,8 @@ void ResolveChimericSequences(
 
 void FindOverlapsAndRepetetiveRegions(
     const std::shared_ptr<thread_pool::ThreadPool>& thread_pool,
-    ram::MinimizerEngine& minimizerEngine, double freq, std::uint8_t kmer_len, double identity,
-    const std::vector<std::unique_ptr<Pile>>& piles,
+    ram::MinimizerEngine& minimizer_engine, double freq, std::uint8_t kmer_len,
+    double identity, const std::vector<std::unique_ptr<Pile>>& piles,
     std::vector<std::vector<biosoup::Overlap>>& overlaps,
     std::vector<std::unique_ptr<biosoup::NucleicAcid>>& sequences) {
   biosoup::Timer timer;
@@ -337,7 +396,6 @@ void FindOverlapsAndRepetetiveRegions(
       sequences_map[sequences[i]->id] = i;
     }
   }
-  
 
   std::uint32_t s = 0;
   for (std::uint32_t i = 0; i < sequences.size(); ++i) {
@@ -359,7 +417,7 @@ void FindOverlapsAndRepetetiveRegions(
 
     timer.Start();
 
-    minimizerEngine.Minimize(sequences.begin() + j, sequences.begin() + i + 1);
+    minimizer_engine.Minimize(sequences.begin() + j, sequences.begin() + i + 1);
 
     std::cerr << "[raven::Graph::Construct] minimized " << j << " - " << i + 1
               << " / " << s << " " << std::fixed << timer.Stop() << "s"
@@ -368,19 +426,19 @@ void FindOverlapsAndRepetetiveRegions(
     timer.Start();
 
     std::vector<std::future<std::vector<biosoup::Overlap>>> thread_futures;
-    minimizerEngine.Filter(freq);
+    minimizer_engine.Filter(freq);
     for (std::uint32_t k = 0; k < i + 1; ++k) {
       thread_futures.emplace_back(thread_pool->Submit(
           [&](std::uint32_t i) -> std::vector<biosoup::Overlap> {
             std::vector<std::uint32_t> filtered;
-            auto dst = minimizerEngine.Map(sequences[i],
+            auto dst = minimizer_engine.Map(sequences[i],
                                            true,   // avoid equal
                                            true,   // avoid symmetric
                                            false,  // minhash
                                            &filtered);
             piles[sequences[i]->id]->AddKmers(filtered, kmer_len,
                                               sequences[i]);  // NOLINT
-          
+
             if (identity != 0) {
               std::uint32_t k = 0;
               for (std::uint32_t j = 0; j < dst.size(); ++j) {
@@ -391,26 +449,25 @@ void FindOverlapsAndRepetetiveRegions(
                 const auto& jt = dst[j];
 
                 auto lhs = sequences[sequences_map[jt.lhs_id]]->InflateData(
-                    jt.lhs_begin,
-                    jt.lhs_end - jt.lhs_begin);
+                    jt.lhs_begin, jt.lhs_end - jt.lhs_begin);
 
                 auto rhs = sequences[sequences_map[jt.rhs_id]]->InflateData(
-                    jt.rhs_begin,
-                    jt.rhs_end - jt.rhs_begin);
+                    jt.rhs_begin, jt.rhs_end - jt.rhs_begin);
                 if (!jt.strand) {
                   biosoup::NucleicAcid rhs_{"", rhs};
                   rhs_.ReverseAndComplement();
                   rhs = rhs_.InflateData();
                 }
 
-                auto result = edlibAlign(
-                    lhs.c_str(), lhs.size(),
-                    rhs.c_str(), rhs.size(),
-                    edlibDefaultAlignConfig());
+                auto result = edlibAlign(lhs.c_str(), lhs.size(), rhs.c_str(),
+                                         rhs.size(), edlibDefaultAlignConfig());
 
-                auto score = result.status == EDLIB_STATUS_OK ?
-                    1. - static_cast<double>(result.editDistance) / std::max(lhs.size(), rhs.size()) :  // NOLINT
-                    0.;
+                auto score =
+                    result.status == EDLIB_STATUS_OK
+                        ? 1. - static_cast<double>(result.editDistance) /
+                                   std::max(lhs.size(), rhs.size())
+                        :  // NOLINT
+                        0.;
 
                 edlibFreeAlignResult(result);
 
@@ -421,7 +478,7 @@ void FindOverlapsAndRepetetiveRegions(
               }
               dst.resize(k);
             }
-            
+
             return dst;
           },
           k));
@@ -662,10 +719,12 @@ void ConstructGraph(
 
   if (graph.stage == -5) {
     FindOverlapsAndCreatePiles(thread_pool, minimizerEngine, sequences,
-                               cfg.freq, graph.piles, overlaps, cfg.kMaxNumOverlaps);
+                               cfg.freq, graph.piles, overlaps,
+                               cfg.kMaxNumOverlaps);
     TrimAndAnnotatePiles(thread_pool, graph.piles, overlaps);
 
-    ResolveContainedReads(graph.piles, overlaps, sequences, thread_pool, cfg.identity);
+    ResolveContainedReads(graph.piles, overlaps, sequences, thread_pool,
+                          cfg.identity);
     ResolveChimericSequences(thread_pool, graph.piles, overlaps, sequences);
 
     ++graph.stage;
@@ -680,8 +739,8 @@ void ConstructGraph(
 
   if (graph.stage == -4) {
     FindOverlapsAndRepetetiveRegions(thread_pool, minimizerEngine, cfg.freq,
-                                     cfg.kmer_len, cfg.identity, graph.piles, overlaps,
-                                     sequences);
+                                     cfg.kmer_len, cfg.identity, graph.piles,
+                                     overlaps, sequences);
     ResolveRepeatInducedOverlaps(thread_pool, graph.piles, overlaps, sequences);
 
     ConstructAssemblyGraph(graph, graph.piles, overlaps, sequences);
diff --git a/RavenLib/src/io.cc b/RavenLib/src/io.cc
index 0ac451c..8845ddc 100644
--- a/RavenLib/src/io.cc
+++ b/RavenLib/src/io.cc
@@ -1,43 +1,99 @@
 #include "raven/io.h"
 
-#include <iostream>
+#include <array>
+#include <fstream>
+#include <functional>
+#include <limits>
+#include <string_view>
+
+#include "bioparser/fasta_parser.hpp"
+#include "bioparser/fastq_parser.hpp"
 
 namespace raven {
 
-std::unique_ptr<bioparser::Parser<biosoup::NucleicAcid>> CreateParser(
-    const std::string& path) {
-  auto is_suffix = [](const std::string_view s, const std::string_view suff) {
-    return s.size() < suff.size()
-               ? false
-               : s.compare(s.size() - suff.size(), suff.size(), suff) == 0;
-  };
-
-  if (is_suffix(path, ".fasta") || is_suffix(path, ".fa") ||
-      is_suffix(path, ".fasta.gz") || is_suffix(path, ".fa.gz")) {
-    try {
+namespace detail {
+
+static constexpr auto kFastaSuffxies =
+    std::array<char const*, 4>{".fasta", "fasta.gz", ".fa", ".fa.gz"};
+
+static constexpr auto kFastqSuffixes =
+    std::array<char const*, 4>{".fastq", ".fastq.gz", ".fq", ".fa.gz"};
+
+static auto IsSuffixFor(std::string_view const suffix,
+                        std::string_view const query) -> bool {
+  return suffix.length() <= query.length()
+             ? suffix == query.substr(query.length() - suffix.length())
+             : false;
+}
+
+static std::unique_ptr<bioparser::Parser<biosoup::NucleicAcid>> CreateParser(
+    const std::filesystem::path& path) {
+  using namespace std::placeholders;
+  if (std::filesystem::exists(path)) {
+    if (std::any_of(kFastaSuffxies.cbegin(), kFastaSuffxies.cend(),
+                    std::bind(IsSuffixFor, _1, path.c_str()))) {
       return bioparser::Parser<biosoup::NucleicAcid>::Create<
-          bioparser::FastaParser>(path);  // NOLINT
-    } catch (const std::invalid_argument& exception) {
-      std::cerr << exception.what() << std::endl;
-      return nullptr;
+          bioparser::FastaParser>(path.c_str());
+    } else if (std::any_of(kFastqSuffixes.cbegin(), kFastqSuffixes.cend(),
+                           std::bind(IsSuffixFor, _1, path.c_str()))) {
+      return bioparser::Parser<biosoup::NucleicAcid>::Create<
+          bioparser::FastqParser>(path.c_str());
     }
   }
-  if (is_suffix(path, ".fastq") || is_suffix(path, ".fq") ||
-      is_suffix(path, ".fastq.gz") || is_suffix(path, ".fq.gz")) {
-    try {
-      return bioparser::Parser<biosoup::NucleicAcid>::Create<
-          bioparser::FastqParser>(path);  // NOLINT
-    } catch (const std::invalid_argument& exception) {
-      std::cerr << exception.what() << std::endl;
-      return nullptr;
+
+  throw std::invalid_argument(
+      "[raven::detail::CreateParser] invalid file path: " + path.string());
+}
+
+}  // namespace detail
+
+std::vector<std::unique_ptr<biosoup::NucleicAcid>> LoadSequences(
+    const std::filesystem::path& path) {
+  auto sparser = detail::CreateParser(path);
+
+  return sparser->Parse(std::numeric_limits<std::uint64_t>::max());
+}
+
+std::vector<std::unique_ptr<biosoup::NucleicAcid>> LoadSequences(
+    std::shared_ptr<thread_pool::ThreadPool> thread_pool,
+    const std::vector<std::filesystem::path>& paths) {
+  auto dst = std::vector<std::unique_ptr<biosoup::NucleicAcid>>();
+
+  {
+    auto parse_futures = std::vector<
+        std::future<std::vector<std::unique_ptr<biosoup::NucleicAcid>>>>();
+
+    auto parse_async = [&thread_pool](const std::filesystem::path& path)
+        -> std::future<std::vector<std::unique_ptr<biosoup::NucleicAcid>>> {
+      return thread_pool->Submit(
+          [](const std::filesystem::path& path)
+              -> std::vector<std::unique_ptr<biosoup::NucleicAcid>> {
+            return LoadSequences(path);
+          },
+          std::cref(path));
+    };
+
+    std::transform(paths.cbegin(), paths.cend(),
+                   std::back_inserter(parse_futures), parse_async);
+
+    for (auto& pf : parse_futures) {
+      auto sequences = pf.get();
+      dst.insert(dst.end(), std::make_move_iterator(sequences.begin()),
+                 std::make_move_iterator(sequences.end()));
     }
+
+    std::sort(dst.begin(), dst.end(),
+              [](const std::unique_ptr<biosoup::NucleicAcid>& lhs,
+                 const std::unique_ptr<biosoup::NucleicAcid>& rhs) -> bool {
+                return lhs->id < rhs->id;
+              });
   }
 
-  std::cerr << "[raven::CreateParser] error: file " << path
-            << " has unsupported format extension (valid extensions: .fasta, "
-            << ".fasta.gz, .fa, .fa.gz, .fastq, .fastq.gz, .fq, .fq.gz)"
-            << std::endl;
-  return nullptr;
+  decltype(dst)(std::make_move_iterator(dst.begin()),
+                std::make_move_iterator(dst.end()))
+      .swap(dst);
+
+  return dst;
 }
 
 }  // namespace raven