diff --git a/cpp/include/seqwin/fasta_reader.hpp b/cpp/include/seqwin/fasta_reader.hpp
index 0a5637d..9db3b3c 100644
--- a/cpp/include/seqwin/fasta_reader.hpp
+++ b/cpp/include/seqwin/fasta_reader.hpp
@@ -6,16 +6,30 @@
 
 namespace seqwin {
 
-struct FastaRecord
-{
-  std::string id;
-  std::string sequence;
+struct FastaRecord {
+    std::string id;
+    std::string sequence;
 };
 
-std::vector<FastaRecord>
-read_fasta(const std::string& assembly_path);
+/**
+ * @brief Read records from a plain-text or gzip-compressed FASTA file.
+ *
+ * @param assembly_path Path to a FASTA or `.gz` FASTA file.
+ * @return Vector of parsed FASTA records.
+ * @throws `std::runtime_error` If the file cannot be opened or the FASTA format is invalid.
+ */
+std::vector<FastaRecord> read_fasta(const std::string& assembly_path);
 
-std::size_t
-est_kmer_number(const std::vector<std::string>& assembly_paths, std::size_t windowsize);
+/**
+ * @brief Estimate the number of minimizers based on the size of assembly files.
+ *
+ * @param assembly_paths Paths to assembly FASTA files (plain or gzipped).
+ * @param windowsize Minimizer window size used for the estimate.
+ * @return Estimated total number of minimizers for all assembly files.
+ */
+std::size_t est_kmer_number(
+    const std::vector<std::string>& assembly_paths,
+    std::size_t windowsize
+);
 
 } // namespace seqwin
diff --git a/cpp/include/seqwin/graph.hpp b/cpp/include/seqwin/graph.hpp
index 50c870d..c5c0934 100644
--- a/cpp/include/seqwin/graph.hpp
+++ b/cpp/include/seqwin/graph.hpp
@@ -9,35 +9,85 @@
 
 namespace seqwin {
 
+/**
+ * @brief Location metadata for a minimizer.
+ */
 struct Kmer {
+    /** 0-based position of the minimizer within its FASTA record. */
     std::uint32_t pos;
+    /** 0-based index of the FASTA record within the assembly. */
     std::uint16_t record_idx;
+    /** Assembly index assigned to the source assembly. */
     std::uint16_t assembly_idx;
 };
 
+/**
+ * @brief Minimizer graph node for one unique minimizer hash.
+ *
+ * The `[start, stop)` range identifies minimizers with this hash.
+ * Before merging, it selects entries in `ThreadGraph.idx`, whose values are indices into `ThreadGraph.kmers`.
+ * After merging, it indexes directly into the parallel `Graph.kmers` and `Graph.idx` arrays.
+ */
 struct Node {
+    /** Hash value of the minimizer represented by this node. */
     std::uint64_t hash;
+    /** Start of the half-open range for this node's minimizer entries. */
     std::uint64_t start;
+    /** End of the half-open range for this node's minimizer entries. */
     std::uint64_t stop;
+    /** Number of target assemblies containing this minimizer hash. */
     std::uint32_t n_tar;
+    /** Number of non-target assemblies containing this minimizer hash. */
     std::uint32_t n_neg;
-    double penalty; // Used as thread_id before merging from different threads
+    /**
+     * Node penalty score used for downstream graph filtering (set to 0.0).
+     * Temporarily stores thread ID before thread-local graphs are merged.
+     */
+    double penalty;
 };
 
+/**
+ * @brief Undirected weighted edge between two minimizers.
+ */
 struct Edge {
+    /** Smaller endpoint hash of the undirected edge. */
     std::uint64_t first;
+    /** Larger endpoint hash of the undirected edge. */
     std::uint64_t second;
+    /** Number of assemblies where the endpoints are adjacent. */
     std::uint64_t weight;
 };
 
+/**
+ * @brief Container for the minimizer graph returned by `build()`.
+ */
 struct Graph {
+    /** Minimizer occurrences in all assemblies, grouped and sorted by hash. */
     NoInitArray<Kmer> kmers;
+    /** Parallel to `kmers` and stores each minimizer's original generation index, ordered by genomic position. */
     NoInitArray<std::uint64_t> idx;
+    /** Sorted by hash. */
     NoInitArray<Node> nodes;
+    /** Sorted by hash. */
     NoInitArray<Edge> edges;
+    /** FASTA record IDs of each assembly. */
     std::vector<std::vector<std::string>> ids_by_assembly;
 };
 
+/**
+ * @brief Build a minimizer graph from assembly FASTA files.
+ *
+ * `assembly_paths`, `assembly_indices`, and `is_targets` are parallel lists.
+ *
+ * @param assembly_paths Paths to input assemblies in FASTA format (plain or gzipped).
+ * @param kmerlen K-mer length for minimizer sketch.
+ * @param windowsize Window size for minimizer sketch.
+ * @param assembly_indices Assembly indices, expected to be `0..N-1`.
+ * @param is_targets Whether each assembly is a target assembly.
+ * @param n_cpu Number of worker threads to use.
+ * @return Minimizer graph.
+ * @throws `std::runtime_error` If input sizes are inconsistent or indices exceed supported ranges.
+ */
 Graph build(
     const std::vector<std::string>& assembly_paths,
     std::size_t kmerlen,
diff --git a/cpp/include/seqwin/helpers.hpp b/cpp/include/seqwin/helpers.hpp
index b881c3f..f5e1993 100644
--- a/cpp/include/seqwin/helpers.hpp
+++ b/cpp/include/seqwin/helpers.hpp
@@ -9,28 +9,54 @@
 
 namespace seqwin {
 
-/**
- * ThreadPool.parallel_for(n_items, fn) splits [0, n_items) into contiguous
- * chunks and calls fn(start, end, worker_id) once per non-empty chunk.
- * Callers should let parallel_for() choose chunk boundaries.
- */
 class ThreadPool;
 
+/**
+ * @brief Partial minimizer graph built by one worker thread.
+ */
 struct ThreadGraph {
-    std::vector<Kmer> kmers; // Needs push_back() support (exact kmer count is not known)
+    /**
+     * Minimizer occurrences generated by this worker, stored in genomic position order.
+     * Uses `std::vector` for `push_back()` support (exact kmer count is not known).
+     */
+    std::vector<Kmer> kmers;
+    /**
+     * Indices into `kmers`, grouped by minimizer hash.
+     * The order is stable within each hash group, preserving the original genomic
+     * position order of minimizers with the same hash.
+     */
     NoInitArray<std::uint64_t> idx;
+    /** Unsorted. */
     NoInitArray<Node> nodes;
+    /** Unsorted. */
     NoInitArray<Edge> edges;
+    /** FASTA record IDs of each assembly. */
     std::vector<std::vector<std::string>> ids_by_assembly;
+    /** Total number of minimizers generated by this worker. */
     std::size_t n_kmers = 0;
+    /** Global index of the first assembly assigned to this worker. */
     std::size_t start_assembly = 0;
 };
 
+/**
+ * @brief Emit a message through Python's logging module.
+ *
+ * @param message Message to log.
+ * @param level Logging level: `debug`, `info`, `warning`, `error`, or `critical`.
+ */
 void log_python(
     const std::string& message,
     const std::string& level = "info"
 );
 
+/**
+ * @brief Merge thread-local minimizer graphs into a single graph.
+ *
+ * @param graphs Thread-local graphs produced during parallel graph construction.
+ * @param n_assemblies Total number of input assemblies.
+ * @param pool Thread pool used for parallel merge steps.
+ * @return Fully merged minimizer graph.
+ */
 Graph merge_thread_graphs(
     std::vector<ThreadGraph>& graphs,
     std::size_t n_assemblies,
diff --git a/cpp/include/seqwin/no_init_array.hpp b/cpp/include/seqwin/no_init_array.hpp
index 5130705..ee30ea1 100644
--- a/cpp/include/seqwin/no_init_array.hpp
+++ b/cpp/include/seqwin/no_init_array.hpp
@@ -7,7 +7,7 @@
 namespace seqwin {
 
 /**
- * A fixed-size owning array.
+ * @brief Fixed-size owning array that avoids value-initializing elements.
  *
  * Unlike `std::vector<T>(n) or std::make_unique<T[]>(n)`, this class allocates
  * with `new T[n]`. For scalar and trivially default-initialized element types,
diff --git a/cpp/include/seqwin/thread_pool.hpp b/cpp/include/seqwin/thread_pool.hpp
index 7e42506..dba4672 100644
--- a/cpp/include/seqwin/thread_pool.hpp
+++ b/cpp/include/seqwin/thread_pool.hpp
@@ -13,6 +13,9 @@
 
 namespace seqwin {
 
+/**
+ * @brief Simple fixed-size worker pool for parallel range processing.
+ */
 class ThreadPool {
 public:
     explicit ThreadPool(std::size_t n_workers)
@@ -45,6 +48,16 @@ class ThreadPool {
 
     std::size_t size() const noexcept { return n_workers_; }
 
+    /**
+     * @brief Run a function over contiguous chunks of `[0, n_items)`.
+     *
+     * The callable receives `(start, end, worker_id)`. Any exception thrown by a
+     * worker is captured and rethrown on the calling thread.
+     *
+     * @tparam Fn Callable type.
+     * @param n_items Number of items in the range.
+     * @param fn Function invoked once for each non-empty chunk.
+     */
     template <typename Fn>
     void parallel_for(std::size_t n_items, Fn&& fn)
     {
diff --git a/cpp/src/bindings/python_bindings.cpp b/cpp/src/bindings/python_bindings.cpp
index db10742..e0ff226 100644
--- a/cpp/src/bindings/python_bindings.cpp
+++ b/cpp/src/bindings/python_bindings.cpp
@@ -71,22 +71,22 @@ PYBIND11_MODULE(_core, m) {
                 owner->edges.data(),
                 capsule
             );
-            py::list all_idx_to_id;
+            py::list ids_by_assembly;
             for (const auto& ids : owner->ids_by_assembly) {
-                py::tuple idx_to_id(ids.size());
+                py::tuple ids_tuple(ids.size());
                 for (std::size_t i = 0; i < ids.size(); ++i) {
-                    idx_to_id[i] = ids[i];
+                    ids_tuple[i] = ids[i];
                 }
-                all_idx_to_id.append(idx_to_id);
+                ids_by_assembly.append(ids_tuple);
             }
 
-            return py::make_tuple(kmers, idx, nodes, edges, all_idx_to_id);
+            return py::make_tuple(kmers, idx, nodes, edges, ids_by_assembly);
         },
         py::arg("assembly_paths"),
         py::arg("kmerlen"),
         py::arg("windowsize"),
-        py::arg("assembly_idx"),
-        py::arg("is_target"),
+        py::arg("assembly_indices"),
+        py::arg("is_targets"),
         py::arg("n_cpu") = 1
     );
 
diff --git a/cpp/src/seqwin/fasta_reader.cpp b/cpp/src/seqwin/fasta_reader.cpp
index c6fb095..720dc32 100644
--- a/cpp/src/seqwin/fasta_reader.cpp
+++ b/cpp/src/seqwin/fasta_reader.cpp
@@ -20,220 +20,219 @@ constexpr std::size_t plain_fasta_seq_len_per_byte = 1;
 constexpr std::size_t gz_fasta_seq_len_per_byte = 4;
 constexpr std::size_t gz_read_buf_size = 1U << 16;
 
-bool
-ends_with(std::string_view text, std::string_view suffix)
+bool ends_with(std::string_view text, std::string_view suffix)
 {
-  return text.size() >= suffix.size() &&
-         text.substr(text.size() - suffix.size()) == suffix;
+    return text.size() >= suffix.size() &&
+           text.substr(text.size() - suffix.size()) == suffix;
 }
 
-std::string
-extract_id(std::string_view header)
+std::string extract_id(std::string_view header)
 {
-  const std::size_t id_end = header.find_first_of(" \t\n\r\f\v");
-  if (id_end == std::string_view::npos) {
-    return std::string(header);
-  }
-  return std::string(header.substr(0, id_end));
+    const std::size_t id_end = header.find_first_of(" \t\n\r\f\v");
+    if (id_end == std::string_view::npos) {
+        return std::string(header);
+    }
+    return std::string(header.substr(0, id_end));
 }
 
-bool
-is_ascii_whitespace_only(std::string_view text)
+bool is_ascii_whitespace_only(std::string_view text)
 {
-  return text.find_first_not_of(" \t\n\r\f\v") == std::string_view::npos;
+    return text.find_first_not_of(" \t\n\r\f\v") == std::string_view::npos;
 }
 
 template<typename NextLine>
-std::vector<FastaRecord>
-read_fasta_core(NextLine&& next_line)
+std::vector<FastaRecord> read_fasta_core(NextLine&& next_line)
 {
-  std::vector<FastaRecord> records;
+    std::vector<FastaRecord> records;
 
-  std::string line;
-  FastaRecord current;
-  bool have_current = false;
+    std::string line;
+    FastaRecord current;
+    bool have_current = false;
 
-  while (next_line(line)) {
-    if (!line.empty() && line.back() == '\r') {
-      line.pop_back();
-    }
+    while (next_line(line)) {
+        if (!line.empty() && line.back() == '\r') {
+            line.pop_back();
+        }
 
-    if (line.empty() || is_ascii_whitespace_only(line)) {
-      continue;
-    }
-    if (line[0] == '>') {
-      if (have_current) {
-        records.push_back(std::move(current));
-        current = FastaRecord{};
-      }
-      current.id = extract_id(line.substr(1));
-      current.sequence.clear();
-      have_current = true;
-      continue;
-    }
+        if (line.empty() || is_ascii_whitespace_only(line)) {
+            continue;
+        }
+        if (line[0] == '>') {
+            if (have_current) {
+                records.push_back(std::move(current));
+                current = FastaRecord{};
+            }
+            current.id = extract_id(line.substr(1));
+            current.sequence.clear();
+            have_current = true;
+            continue;
+        }
 
-    if (!have_current) {
-      throw std::runtime_error("Invalid FASTA: sequence encountered before header");
-    }
+        if (!have_current) {
+            throw std::runtime_error("Invalid FASTA: sequence encountered before header");
+        }
 
-    const auto first_ws = line.find_first_of(" \t\n\r\f\v");
-    if (first_ws == std::string::npos) {
-      current.sequence.append(line);
-      continue;
-    }
+        const auto first_ws = line.find_first_of(" \t\n\r\f\v");
+        if (first_ws == std::string::npos) {
+            current.sequence.append(line);
+            continue;
+        }
 
-    const std::size_t needed = current.sequence.size() + line.size();
-    if (needed > current.sequence.capacity()) {
-      const std::size_t grown = std::max(current.sequence.capacity() * 2, needed);
-      current.sequence.reserve(grown);
-    }
-    for (char c : line) {
-      if (!std::isspace(static_cast<unsigned char>(c))) {
-        current.sequence.push_back(c);
-      }
+        const std::size_t needed = current.sequence.size() + line.size();
+        if (needed > current.sequence.capacity()) {
+            const std::size_t grown = std::max(current.sequence.capacity() * 2, needed);
+            current.sequence.reserve(grown);
+        }
+        for (char c : line) {
+            if (!std::isspace(static_cast<unsigned char>(c))) {
+                current.sequence.push_back(c);
+            }
+        }
     }
-  }
 
-  if (have_current) {
-    records.push_back(std::move(current));
-  }
+    if (have_current) {
+        records.push_back(std::move(current));
+    }
 
-  return records;
+    return records;
 }
 
-std::vector<FastaRecord>
-read_plain_fasta(const std::string& assembly_path)
+std::vector<FastaRecord> read_plain_fasta(const std::string& assembly_path)
 {
-  std::ifstream in(assembly_path);
-  if (!in) {
-    throw std::runtime_error("Unable to open FASTA: " + assembly_path);
-  }
-
-  return read_fasta_core([&in](std::string& line) -> bool {
-    return static_cast<bool>(std::getline(in, line));
-  });
+    std::ifstream in(assembly_path);
+    if (!in) {
+        throw std::runtime_error("Unable to open FASTA: " + assembly_path);
+    }
+
+    return read_fasta_core([&in](std::string& line) -> bool {
+        return static_cast<bool>(std::getline(in, line));
+    });
 }
 
-std::vector<FastaRecord>
-read_gz_fasta(const std::string& assembly_path)
+std::vector<FastaRecord> read_gz_fasta(const std::string& assembly_path)
 {
-  gzFile raw_gz = gzopen(assembly_path.c_str(), "rb");
-  if (raw_gz == nullptr) {
-    throw std::runtime_error("Unable to open gzip FASTA: " + assembly_path);
-  }
-
-  struct GzCloser
-  {
-    void operator()(gzFile_s* f) const
-    {
-      if (f != nullptr) {
-        gzclose(reinterpret_cast<gzFile>(f));
-      }
+    gzFile raw_gz = gzopen(assembly_path.c_str(), "rb");
+    if (raw_gz == nullptr) {
+        throw std::runtime_error("Unable to open gzip FASTA: " + assembly_path);
     }
-  };
 
-  std::unique_ptr<gzFile_s, GzCloser> gz(reinterpret_cast<gzFile_s*>(raw_gz));
+    struct GzCloser {
+        void operator()(gzFile_s* f) const
+        {
+            if (f != nullptr) {
+                gzclose(reinterpret_cast<gzFile>(f));
+            }
+        }
+    };
 
-  std::string carry;
-  std::size_t line_start = 0;
-  std::size_t search_start = 0;
-  bool at_eof = false;
+    std::unique_ptr<gzFile_s, GzCloser> gz(reinterpret_cast<gzFile_s*>(raw_gz));
 
-  auto fill_buffer = [&]() {
-    if (at_eof) {
-      return;
-    }
+    std::string carry;
+    std::size_t line_start = 0;
+    std::size_t search_start = 0;
+    bool at_eof = false;
 
-    std::array<char, gz_read_buf_size> buf{};
-    int bytes =
-      gzread(reinterpret_cast<gzFile>(gz.get()), buf.data(), static_cast<unsigned int>(buf.size()));
-    if (bytes < 0) {
-      int errnum = 0;
-      const char* err = gzerror(reinterpret_cast<gzFile>(gz.get()), &errnum);
-      throw std::runtime_error(std::string("gzip read error: ") + (err ? err : "unknown"));
-    }
-    if (bytes == 0) {
-      at_eof = true;
-      return;
-    }
-    carry.append(buf.data(), static_cast<std::size_t>(bytes));
-  };
-
-  auto records = read_fasta_core([&](std::string& line) -> bool {
-    for (;;) {
-      const std::size_t nl = carry.find('\n', search_start);
-      if (nl != std::string::npos) {
-        std::size_t end = nl;
-        if (end > line_start && carry[end - 1] == '\r') {
-          --end;
+    auto fill_buffer = [&]() {
+        if (at_eof) {
+            return;
+        }
+
+        std::array<char, gz_read_buf_size> buf{};
+        int bytes = gzread(
+            reinterpret_cast<gzFile>(gz.get()),
+            buf.data(),
+            static_cast<unsigned int>(buf.size())
+        );
+        if (bytes < 0) {
+            int errnum = 0;
+            const char* err = gzerror(reinterpret_cast<gzFile>(gz.get()), &errnum);
+            throw std::runtime_error(
+                std::string("gzip read error: ") + (err ? err : "unknown")
+            );
         }
-        line.assign(carry.data() + line_start, end - line_start);
-        line_start = nl + 1;
-        search_start = line_start;
-
-        if (line_start > gz_read_buf_size && line_start >= carry.size() / 2) {
-          carry.erase(0, line_start);
-          search_start -= line_start;
-          line_start = 0;
+        if (bytes == 0) {
+            at_eof = true;
+            return;
         }
-        return true;
-      }
-
-      search_start = carry.size();
-
-      if (at_eof) {
-        if (line_start < carry.size()) {
-          std::size_t end = carry.size();
-          if (end > line_start && carry[end - 1] == '\r') {
-            --end;
-          }
-          line.assign(carry.data() + line_start, end - line_start);
-          line_start = carry.size();
-          search_start = line_start;
-          return true;
+        carry.append(buf.data(), static_cast<std::size_t>(bytes));
+    };
+
+    auto records = read_fasta_core([&](std::string& line) -> bool {
+        for (;;) {
+            const std::size_t nl = carry.find('\n', search_start);
+            if (nl != std::string::npos) {
+                std::size_t end = nl;
+                if (end > line_start && carry[end - 1] == '\r') {
+                    --end;
+                }
+                line.assign(carry.data() + line_start, end - line_start);
+                line_start = nl + 1;
+                search_start = line_start;
+
+                if (line_start > gz_read_buf_size && line_start >= carry.size() / 2) {
+                    carry.erase(0, line_start);
+                    search_start -= line_start;
+                    line_start = 0;
+                }
+                return true;
+            }
+
+            search_start = carry.size();
+
+            if (at_eof) {
+                if (line_start < carry.size()) {
+                    std::size_t end = carry.size();
+                    if (end > line_start && carry[end - 1] == '\r') {
+                        --end;
+                    }
+                    line.assign(carry.data() + line_start, end - line_start);
+                    line_start = carry.size();
+                    search_start = line_start;
+                    return true;
+                }
+                return false;
+            }
+
+            fill_buffer();
+            if (line_start > gz_read_buf_size && line_start >= carry.size() / 2) {
+                carry.erase(0, line_start);
+                search_start -= line_start;
+                line_start = 0;
+            }
         }
-        return false;
-      }
-
-      fill_buffer();
-      if (line_start > gz_read_buf_size && line_start >= carry.size() / 2) {
-        carry.erase(0, line_start);
-        search_start -= line_start;
-        line_start = 0;
-      }
-    }
-  });
+    });
 
-  return records;
+    return records;
 }
 
 } // namespace
 
-std::vector<FastaRecord>
-read_fasta(const std::string& assembly_path)
+std::vector<FastaRecord> read_fasta(const std::string& assembly_path)
 {
-  if (ends_with(assembly_path, ".gz")) {
-    return read_gz_fasta(assembly_path);
-  }
-  return read_plain_fasta(assembly_path);
+    if (ends_with(assembly_path, ".gz")) {
+        return read_gz_fasta(assembly_path);
+    }
+    return read_plain_fasta(assembly_path);
 }
 
-std::size_t
-est_kmer_number(const std::vector<std::string>& assembly_paths, std::size_t windowsize)
-{
-  // Reserve-only heuristic, not correctness-critical.
-  std::size_t est_total_seq_len = 0;
-  for (const auto& assembly_path : assembly_paths) {
-    const std::size_t seq_len_per_byte =
-      ends_with(assembly_path, ".gz") ? gz_fasta_seq_len_per_byte : plain_fasta_seq_len_per_byte;
-    std::error_code ec;
-    const auto file_bytes = std::filesystem::file_size(assembly_path, ec);
-    if (ec) {
-      continue;
+std::size_t est_kmer_number(
+    const std::vector<std::string>& assembly_paths,
+    std::size_t windowsize
+) {
+    // Reserve-only heuristic, not correctness-critical.
+    std::size_t est_total_seq_len = 0;
+    for (const auto& assembly_path : assembly_paths) {
+        const std::size_t seq_len_per_byte = ends_with(assembly_path, ".gz")
+            ? gz_fasta_seq_len_per_byte
+            : plain_fasta_seq_len_per_byte;
+        std::error_code ec;
+        const auto file_bytes = std::filesystem::file_size(assembly_path, ec);
+        if (ec) {
+            continue;
+        }
+        est_total_seq_len += static_cast<std::size_t>(file_bytes * seq_len_per_byte);
     }
-    est_total_seq_len += static_cast<std::size_t>(file_bytes * seq_len_per_byte);
-  }
-  return (2 * est_total_seq_len) / (windowsize + 1);
+    return (2 * est_total_seq_len) / (windowsize + 1);
 }
 
 } // namespace seqwin
diff --git a/cpp/src/seqwin/graph.cpp b/cpp/src/seqwin/graph.cpp
index 881a54b..c4f87a1 100644
--- a/cpp/src/seqwin/graph.cpp
+++ b/cpp/src/seqwin/graph.cpp
@@ -53,6 +53,7 @@ ThreadGraph build_worker(
     std::size_t end_assembly,
     std::size_t thread_id
 ) {
+    // Estimate total minimizer count in all assemblies
     const auto n_kmers_est = seqwin::est_kmer_number(
         std::vector<std::string>(
             assembly_paths.begin() + static_cast<std::ptrdiff_t>(start_assembly),
@@ -66,7 +67,7 @@ ThreadGraph build_worker(
     graph.ids_by_assembly.resize(end_assembly - start_assembly);
     graph.start_assembly = start_assembly;
 
-    std::vector<std::uint64_t> hashes;
+    std::vector<std::uint64_t> hashes; // Used to build ThreadGraph.idx
     hashes.reserve(n_kmers_est);
 
     // Reserving for unordered_map will actually allocate physical memory
@@ -78,38 +79,34 @@ ThreadGraph build_worker(
 
     for (std::size_t assembly_i = start_assembly; assembly_i < end_assembly; ++assembly_i) {
         const auto assembly_idx = assembly_indices[assembly_i];
-        if (assembly_idx > std::numeric_limits<std::uint16_t>::max()) {
-            throw std::runtime_error("assembly_idx must fit in uint16");
-        }
-        const auto assembly_idx16 = static_cast<std::uint16_t>(assembly_idx);
         const bool is_target = is_targets[assembly_i];
 
         auto records = seqwin::read_fasta(assembly_paths[assembly_i]);
-        auto& record_IDs = graph.ids_by_assembly[assembly_i - start_assembly];
-        record_IDs.resize(records.size());
+        auto& record_ids = graph.ids_by_assembly[assembly_i - start_assembly];
+        record_ids.resize(records.size());
 
         for (std::size_t record_idx = 0; record_idx < records.size(); ++record_idx) {
             if (record_idx > std::numeric_limits<std::uint16_t>::max()) {
-                throw std::runtime_error("record_idx must fit in uint16");
+                throw std::runtime_error("Number of FASTA records exceeds uint16 range");
             }
-            const auto record_idx16 = static_cast<std::uint16_t>(record_idx);
-
             auto& record = records[record_idx];
-            record_IDs[record_idx] = std::move(record.id);
+            record_ids[record_idx] = std::move(record.id);
 
+            // Generate minimizers for the current record
             const auto mins = btllib::minimize_sequence(record.sequence, kmerlen, windowsize);
 
             for (const auto& m : mins) {
                 if (m.pos > std::numeric_limits<std::uint32_t>::max()) {
-                    throw std::runtime_error("minimizer position exceeds uint32 range");
+                    throw std::runtime_error("Minimizer position exceeds uint32 range");
                 }
                 graph.kmers.push_back(Kmer{
                     static_cast<std::uint32_t>(m.pos),
-                    record_idx16,
-                    assembly_idx16
+                    static_cast<std::uint16_t>(record_idx),
+                    static_cast<std::uint16_t>(assembly_idx)
                 });
                 hashes.push_back(m.out_hash);
 
+                // Add this minimizer to an existing node, or create a new node
                 auto [node_it, node_inserted] = node_map.try_emplace(m.out_hash);
                 ++(node_it->second.count);
                 if (node_inserted || node_it->second.last_seen_assembly != assembly_i) {
@@ -123,10 +120,12 @@ ThreadGraph build_worker(
                 ++graph.n_kmers;
             }
 
+            // Current record is too short for an edge
             if (mins.size() < 2) {
                 continue;
             }
 
+            // Add undirected edges
             for (std::size_t i = 0; i + 1 < mins.size(); ++i) {
                 auto u = mins[i].out_hash;
                 auto v = mins[i + 1].out_hash;
@@ -143,7 +142,7 @@ ThreadGraph build_worker(
         }
     }
 
-    // Create edges first to reduce peak memory
+    // Materialize edges first to reduce peak memory
     graph.edges = NoInitArray<Edge>(edge_map.size());
     std::size_t edge_i = 0;
     for (const auto& [key, state] : edge_map) {
@@ -151,6 +150,7 @@ ThreadGraph build_worker(
     }
     std::unordered_map<EdgeKey, EdgeState, EdgeKeyHash>().swap(edge_map);
 
+    // Build ThreadGraph.idx (grouped by minimizer hash)
     graph.idx = NoInitArray<std::uint64_t>(graph.n_kmers);
     std::uint64_t cursor = 0;
     for (auto& [hash, state] : node_map) {
@@ -195,11 +195,11 @@ Graph build(
     if (assembly_paths.size() != assembly_indices.size() ||
         assembly_paths.size() != is_targets.size()) {
         throw std::runtime_error(
-            "assembly_paths, assembly_idx, and is_target must have the same length");
+            "assembly_paths, assembly_indices, and is_targets must have the same length");
     }
     for (std::size_t i = 0; i < assembly_indices.size(); ++i) {
-        if (assembly_indices[i] != i) {
-            throw std::runtime_error("assembly_indices must be strictly incremental from 0 to N-1");
+        if (assembly_indices[i] > std::numeric_limits<std::uint16_t>::max()) {
+            throw std::runtime_error("Assembly index exceeds uint16 range");
         }
     }
 
@@ -209,7 +209,7 @@ Graph build(
         n_workers = std::min(n_workers, n_assemblies);
     }
 
-    ThreadPool pool(n_workers);
+    ThreadPool pool(n_workers); // Avoid spawning threads every time
     std::vector<ThreadGraph> graphs(n_workers);
 
     const std::size_t base = n_assemblies / n_workers;
diff --git a/cpp/src/seqwin/helpers.cpp b/cpp/src/seqwin/helpers.cpp
index 48799ff..190a94b 100644
--- a/cpp/src/seqwin/helpers.cpp
+++ b/cpp/src/seqwin/helpers.cpp
@@ -13,6 +13,9 @@ namespace py = pybind11;
 namespace seqwin {
 namespace {
 
+/**
+ * @brief Describes a contiguous local `idx` segment and its output position.
+ */
 struct IdxSegment {
     std::size_t thread_id;
     std::size_t local_start;
@@ -64,9 +67,15 @@ NoInitArray<Edge> concat_edges(std::vector<ThreadGraph>& graphs, ThreadPool& poo
 }
 
 /**
- * 1. Parallel LSD radix sort based on hash (stable)
- * 2. Merge nodes with the same hash
- * 3. Build segment metadata for final materialization
+ * @brief Sort and merge nodes with identical hashes.
+ *
+ * 1. Parallel LSD radix sort on node hash.
+ * 2. Merge nodes with the same hash.
+ * 3. Return local `idx` segment metadata to materialize `kmers`.
+ *
+ * @param nodes Node array to sort and merge.
+ * @param pool Thread pool used for radix sorting.
+ * @return Segment metadata for rebuilding `kmers` and `idx`.
  */
 static std::vector<IdxSegment> merge_nodes(
     NoInitArray<Node>& nodes,
diff --git a/src/seqwin/__init__.py b/src/seqwin/__init__.py
index 180a808..652cbf0 100644
--- a/src/seqwin/__init__.py
+++ b/src/seqwin/__init__.py
@@ -3,10 +3,10 @@
 ======
 
 Ultrafast identification of signature sequences in microbial genomes
-- Output sensitive and specific genomic signatures. 
-- High input scalability. 
-- Robust to sequence variations and low-quality assemblies. 
-- Use raw sequences as input. 
+- Output sensitive and specific genomic signatures.
+- High input scalability.
+- Robust to sequence variations and low-quality assemblies.
+- Use raw sequences as input.
 
 Usage:
 ------
diff --git a/src/seqwin/assemblies.py b/src/seqwin/assemblies.py
index 5e8a9d0..f5b0962 100644
--- a/src/seqwin/assemblies.py
+++ b/src/seqwin/assemblies.py
@@ -2,7 +2,7 @@
 Assemblies
 ==========
 
-Create an instance for all input genome assemblies. 
+Create an instance for all input genome assemblies.
 
 Dependencies:
 -------------
@@ -32,7 +32,6 @@
 from io import BufferedWriter
 from time import time
 from queue import Empty
-from collections.abc import Callable
 
 logger = logging.getLogger(__name__)
 
@@ -47,52 +46,52 @@
 from .config import Config, RunState, WORKINGDIR, BLASTCONFIG
 
 _FASTA_EXT = (
-    '.fna', '.fasta', '.fna.gz', '.fasta.gz', 
+    '.fna', '.fasta', '.fna.gz', '.fasta.gz',
     '.fa', '.fas', '.fa.gz', '.fas.gz'
 )
 
 
 class Assemblies(pd.DataFrame):
-    """Package all input genome assemblies as a pandas DataFrame. 
+    """Package all input genome assemblies as a pandas DataFrame.
 
     Attributes:
-        path (pd.Series[Path]): Assembly paths. 
-        is_target (pd.Series[bool]): True for target assemblies. 
-        record_ids (pd.Series[tuple[str, ...] | None]): Record IDs of each assembly. 
+        path (pd.Series[Path]): Assembly paths.
+        is_target (pd.Series[bool]): True for target assemblies.
+        record_ids (pd.Series[tuple[str, ...] | None]): Record IDs of each assembly.
     """
     def __init__(self, tar_paths: list[Path], neg_paths: list[Path]) -> None:
-        """Package all input genome assemblies as a pandas DataFrame. 
+        """Package all input genome assemblies as a pandas DataFrame.
 
         Args:
-            tar_paths (list[Path]): A list of paths to target assemblies. 
-            neg_paths (list[Path]): A list of paths to non-target assemblies. 
+            tar_paths (list[Path]): A list of paths to target assemblies.
+            neg_paths (list[Path]): A list of paths to non-target assemblies.
         """
         data = dict(
-            path=tar_paths + neg_paths, 
-            is_target=[True]*len(tar_paths) + [False]*len(neg_paths), 
+            path=tar_paths + neg_paths,
+            is_target=[True]*len(tar_paths) + [False]*len(neg_paths),
             record_ids=None
         )
         super().__init__(data)
 
     def mash(self, kmerlen: int, sketchsize: int, out_path: Path, overwrite: bool, n_cpu: int) -> NDArray[np.floating]:
-        """Calculate the Jaccard indices of all assembly pairs with Mash. 
+        """Calculate the Jaccard indices of all assembly pairs with Mash.
 
         Args:
-            kmerlen (int): K-mer length for `mash sketch`. 
-            sketchsize (int): Sketch size for `mash sketch`. 
-            out_path (Path): Output path for the Mash sketch file (.msh). 
-            overwrite (bool): If True, overwrite the existing Mash sketch file. 
-            n_cpu (int): Number of processes to run in parallel. 
+            kmerlen (int): K-mer length for `mash sketch`.
+            sketchsize (int): Sketch size for `mash sketch`.
+            out_path (Path): Output path for the Mash sketch file (.msh).
+            overwrite (bool): If True, overwrite the existing Mash sketch file.
+            n_cpu (int): Number of processes to run in parallel.
 
         Returns:
-            NDArray[np.floating]: A matrix of Jaccard indices of all assembly pairs. 
+            NDArray[np.floating]: A matrix of Jaccard indices of all assembly pairs.
         """
         mash_sketch = sketch(
-            self.path.tolist(), 
-            kmerlen=kmerlen, 
-            sketchsize=sketchsize, 
-            out_path=out_path, 
-            overwrite=overwrite, 
+            self.path.tolist(),
+            kmerlen=kmerlen,
+            sketchsize=sketchsize,
+            out_path=out_path,
+            overwrite=overwrite,
             n_cpu=n_cpu
         )
         return np.array(
@@ -100,21 +99,21 @@ def mash(self, kmerlen: int, sketchsize: int, out_path: Path, overwrite: bool, n
         ).reshape(len(self), len(self))
 
     def fetch_seq(self, loc: pd.DataFrame, n_cpu: int) -> pd.Series:
-        """Fetch the actual sequences for a DataFrame of assembly locations. 
-        - Fetching the sequence of each location one by one is slow, since it needs layers of indices to 
-        access the actual sequence (assembly, record, start and stop). 
-        - To solve this, rows from the same assembly are grouped together, 
-        and different groups are fetched in parallel. 
+        """Fetch the actual sequences for a DataFrame of assembly locations.
+        - Fetching the sequence of each location one by one is slow, since it needs layers of indices to
+        access the actual sequence (assembly, record, start and stop).
+        - To solve this, rows from the same assembly are grouped together,
+        and different groups are fetched in parallel.
 
         Args:
-            loc (pd.DataFrame): Assembly locations. Row indices are kept in the returned Series, but the 
-                ordering might be different. To make sure the returned Series has the same order as `loc`, 
-                row indices should be sorted with `ascending=True`. 
-                Required columns: ['assembly_idx', 'record_idx', 'start', 'stop']. 
-            n_cpu (int): Number of processes to run in parallel. 
+            loc (pd.DataFrame): Assembly locations. Row indices are kept in the returned Series, but the
+                ordering might be different. To make sure the returned Series has the same order as `loc`,
+                row indices should be sorted with `ascending=True`.
+                Required columns: ['assembly_idx', 'record_idx', 'start', 'stop'].
+            n_cpu (int): Number of processes to run in parallel.
 
         Returns:
-            pd.Series: A sequence is fetched for each row in `loc`. indices are sorted with `ascending=True`. 
+            pd.Series: A sequence is fetched for each row in `loc`. indices are sorted with `ascending=True`.
         """
         # group sequences by assembly_idx
         loc: dict[int, pd.DataFrame] = dict(tuple(
@@ -130,11 +129,11 @@ def fetch_seq(self, loc: pd.DataFrame, n_cpu: int) -> pd.Series:
 
         # fetch the actual sequences by start and stop in the source sequences
         fetch_seq_args = zip(
-            loc.values(), 
+            loc.values(),
             all_src_paths
         )
         all_seq: pd.Series = pd.concat(
-            mp_wrapper(_fetch_seq, fetch_seq_args, n_cpu, n_jobs=len(loc)), 
+            mp_wrapper(_fetch_seq, fetch_seq_args, n_cpu, n_jobs=len(loc)),
             axis=0
         )
         # sort the returned sequences by the original ordering (before groupby)
@@ -142,21 +141,21 @@ def fetch_seq(self, loc: pd.DataFrame, n_cpu: int) -> pd.Series:
         return all_seq
 
     def makeblastdb(self, prefix: Path, neg_only: bool, overwrite: bool, n_cpu: int) -> Path:
-        """Create a BLAST database for all (or non-target) assemblies. Use native Python streaming and multiprocessing. 
-        - Note: macOS (x64 or ARM) has a hard-wired pipe buffer size of 64kB (vs. 1MB on Linux), so `makeblastdb` will 
+        """Create a BLAST database for all (or non-target) assemblies. Use native Python streaming and multiprocessing.
+        - Note: macOS (x64 or ARM) has a hard-wired pipe buffer size of 64kB (vs. 1MB on Linux), so `makeblastdb` will
         be a lot slower on a Mac when the input is streamed to `stdin`. While on Linux the difference is negligible
-        due to the larger buffer size. 
+        due to the larger buffer size.
 
         Args:
-            prefix (Path): Output directory of the BLAST database. 
-            neg_only (bool): If True, create the BLAST database on non-target assemblies only. 
-            overwrite (bool): If True, overwrite prefix if it already exists. 
-            n_cpu (int): Number of processes to run in parallel. 
+            prefix (Path): Output directory of the BLAST database.
+            neg_only (bool): If True, create the BLAST database on non-target assemblies only.
+            overwrite (bool): If True, overwrite prefix if it already exists.
+            n_cpu (int): Number of processes to run in parallel.
 
         Returns:
-            Path: Path to the BLAST database. 
+            Path: Path to the BLAST database.
         """
-        # NOTE: when the size of the blastdb changes, the evalue of a specific hit also changes. 
+        # NOTE: when the size of the blastdb changes, the evalue of a specific hit also changes.
         # Since the evalue threshold for a blast task is set, this hit might not be included when the blastdb gets larger
         if neg_only:
             logger.info('Creating a BLAST database of non-target assemblies (less sensitive but faster)...')
@@ -175,19 +174,19 @@ def makeblastdb(self, prefix: Path, neg_only: bool, overwrite: bool, n_cpu: int)
         # load fasta files and stream to stdin to save memory
         with mp.Manager() as manager:
             # only Manager().Queue() can be shared between processes (e.g, when used with Pool())
-            # so that different processes can put items in the same queue. 
+            # so that different processes can put items in the same queue.
             queue = manager.Queue(maxsize=BLASTCONFIG.queue_size+n_cpu) # set queue size to limit memory usage
             queue_idx = range(len(df)) # queue index must start from 0 (for df.index, this is not True when neg_only is True)
 
             # create a process for makeblastdb
             makeblastdb_args = [
-                'makeblastdb', 
-                '-title', title, 
-                '-dbtype', 'nucl', 
+                'makeblastdb',
+                '-title', title,
+                '-dbtype', 'nucl',
                 '-out', blastdb
             ]
             proc = subprocess.Popen(
-                makeblastdb_args, stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE, 
+                makeblastdb_args, stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,
                 text=False # use bytes
             )
 
@@ -208,8 +207,8 @@ def makeblastdb(self, prefix: Path, neg_only: bool, overwrite: bool, n_cpu: int)
         # save command, stdout and stderr
         blast_log = prefix / WORKINGDIR.blast_log
         blast_log.write_text('\n'.join((
-            str(makeblastdb_args), 
-            stdout, 
+            str(makeblastdb_args),
+            stdout,
             stderr
         )))
         if proc.returncode != 0:
@@ -218,95 +217,23 @@ def makeblastdb(self, prefix: Path, neg_only: bool, overwrite: bool, n_cpu: int)
         logger.info(f' - BLAST database created: {blastdb}')
         print_time_delta(time()-tik)
         return blastdb
-    
-    def __makeblastdb_cmd(self, title: str, prefix: Path, neg_only: bool, overwrite: bool, n_cpu: int) -> Path:
-        """Create a BLAST database for all or non-target assemblies. Use streaming with subprocess.run(shell=True) to save memory. 
-        `sd` is needed for modifying FASTA content, and `parallel` is needed for multiprocessing. 
-
-        Dependencies:
-        - sd (https://github.com/chmln/sd, https://anaconda.org/conda-forge/sd)
-        - parallel (https://www.gnu.org/software/parallel/, https://anaconda.org/conda-forge/parallel)
-
-        Args:
-            title (str): Name of the BLAST database. 
-            prefix (Path): Output directory of the BLAST database. 
-            neg_only (bool): If True, create the BLAST database on non-target assemblies only. 
-            overwrite (bool): If True, overwrite prefix if it already exists. 
-            n_cpu (int): Number of processes to run in parallel. 
-
-        Returns:
-            Path: Path to the BLAST database. 
-        """
-        # NOTE: when the size of the blastdb changes, the evalue of a specific hit also changes. 
-        # Since the evalue threshold for a blast task is set, this hit might not be included when the blastdb gets larger
-        if neg_only:
-            logger.info('Creating a BLAST database of non-target assemblies (less sensitive but faster)...')
-            df = self[self.is_target == False]
-            title += BLASTCONFIG.neg_only
-        else:
-            logger.info('Creating a BLAST database of all assemblies (more sensitive but slower)...')
-            df = self
-        tik = time()
-
-        # create a folder for BLAST
-        mkdir(prefix, overwrite)
-        blastdb = prefix / title
-
-        # make a tsv of assembly path, assembly index and is_target (input for GNU parallel)
-        assemblies_tsv = '\n'.join(
-            f'{p}\t{i}\t{BLASTCONFIG.bool2str[b]}' 
-            for p, i, b in zip(df.path, df.index, df.is_target)
-        )
-        
-        # GNU parallel command to modify assembly fasta headers (>record_id to >idx@is_target@record_id)
-        # {1}=assembly path, {2}=assembly index, {3}=is_target
-        # use zcat if the fasta file is gzipped, else use cat; then modify the headers with sd
-        # since this command uses shell syntax ([[ ... ]], &&, ||, and |), shell=True is necessary in subprocess
-        modify_fasta = '[[ {1} == *.gz ]] && zcat "{1}" || cat "{1}" | sd "^>" ">{2}%s{3}%s"' \
-            % (BLASTCONFIG.header_sep, BLASTCONFIG.header_sep)
-        
-        # command to create the blast db
-        makeblastdb = f'makeblastdb -title {title} -dbtype nucl -out {blastdb}'
-        
-        # load assembly files in paralle, and use pipes to save memory (shell=True)
-        # tried to use native python pipes with shell=False but failed
-        proc = subprocess.run(
-            rf"parallel -j {n_cpu} --colsep '\t' '{modify_fasta}' | {makeblastdb}", 
-            input=assemblies_tsv, 
-            shell=True, capture_output=True, text=True
-        )
-
-        # save the full command, stdout and stderr
-        blast_log = prefix / WORKINGDIR.blast_log
-        blast_log.write_text('\n'.join((
-            proc.args, 
-            proc.stdout, 
-            proc.stderr, 
-            assemblies_tsv
-        )))
-        if proc.returncode != 0:
-            log_and_raise(RuntimeError, msg=f'Failed to create the BLAST database. For details, please check {blast_log}')
-        
-        logger.info(f' - BLAST database created: {blastdb}')
-        print_time_delta(time()-tik)
-        return blastdb
 
 
 def _add_fasta_to_queue(path: Path, assembly_idx: int, is_target: bool, queue_idx: int, queue: mp.Queue) -> None:
     """
-    1. Add assembly index and is_target into the header lines of an assembly FASTA file. 
-    2. Put the modified FASTA file content, as well as its queue_idx in a queue. 
+    1. Add assembly index and is_target into the header lines of an assembly FASTA file.
+    2. Put the modified FASTA file content, as well as its queue_idx in a queue.
 
-    The use of queue_idx is to make sure that items in the queue can be fetched in order (in `_stream_to_stdin()`), 
-    so that the behavior is exactly the same as reading all FASTA files into memory using a single thread. 
+    The use of queue_idx is to make sure that items in the queue can be fetched in order (in `_stream_to_stdin()`),
+    so that the behavior is exactly the same as reading all FASTA files into memory using a single thread.
 
     Args:
-        path (Path): Path to the assembly FASTA file. 
-        assembly_idx (int): Assembly index. 
-        is_target (bool): True for target assemblies. 
-        queue_idx (int): Used to determine the order of assemblies in the queue. Must start from 0. 
-        queue (mp.Queue): queue_idx and modified FASTA file content are put into this queue. 
-            Use `mp.Manager().Queue()` to share this queue across different processes (e.g., when this function is called by `mp.Pool()`). 
+        path (Path): Path to the assembly FASTA file.
+        assembly_idx (int): Assembly index.
+        is_target (bool): True for target assemblies.
+        queue_idx (int): Used to determine the order of assemblies in the queue. Must start from 0.
+        queue (mp.Queue): queue_idx and modified FASTA file content are put into this queue.
+            Use `mp.Manager().Queue()` to share this queue across different processes (e.g., when this function is called by `mp.Pool()`).
     """
     # read file content as bytes
     if path.suffix == GZIP_EXT:
@@ -326,12 +253,12 @@ def _add_fasta_to_queue(path: Path, assembly_idx: int, is_target: bool, queue_id
 
 
 def _stream_to_stdin(queue: mp.Queue, n_items: int, proc_stdin: BufferedWriter) -> None:
-    """Get items from an indexed queue, and write them to the stdin of a process by the order of their indices. 
+    """Get items from an indexed queue, and write them to the stdin of a process by the order of their indices.
 
     Args:
-        queue (mp.Queue): Each queue item should be a tuple of (idx, data). 
-        n_items (int): Total number of items in the queue. 
-        proc_stdin (io.BufferedWriter): Standard input of a process (e.g., `subprocess.Popen().stdin`). 
+        queue (mp.Queue): Each queue item should be a tuple of (idx, data).
+        n_items (int): Total number of items in the queue.
+        proc_stdin (io.BufferedWriter): Standard input of a process (e.g., `subprocess.Popen().stdin`).
     """
     next_idx = 0
     buffer: dict[int, bytes] = dict()
@@ -352,77 +279,48 @@ def _stream_to_stdin(queue: mp.Queue, n_items: int, proc_stdin: BufferedWriter)
     proc_stdin.flush() # empty stdin but don't close the process
 
 
-def _load_seq(
-    paths: list[Path], 
-    parser: Callable[[Path], tuple[dict[str, str], int]] = load_fasta, 
-    n_cpu: int = 1
-) -> list[dict[str, str]]:
-    """Deprecated. Load assembly sequences from files. 
-
-    Args:
-        paths (list[Path]): A list of paths to the assembly files. 
-        parser (Callable, optional): The function to parse the assembly files. It should return a dict of record id -> 
-            sequence (upper case), and the total number of bases in the assembly. Supported functions: `parsers.load_fasta()`, 
-            `parsers.load_genbank()`. [load_fasta]
-        n_cpu (int, optional): Number of processes to run in parallel. [1]
-
-    Returns:
-        list[dict[str,str]]: Each dict has record IDs as keys and record sequences as values. 
-    """
-    # load assembly files in paralelle
-    n_assemblies = len(paths)
-    logger.info(f' - Loading {n_assemblies} assembly files...')
-    all_seq, all_len = mp_wrapper(
-        parser, paths, n_cpu, 
-        starmap=False, unpack_output=True
-    )
-    total_len = sum(all_len)
-    logger.info(f' - Average assembly size: {total_len/n_assemblies:.0f} bp; total: {total_len} bp')
-    return all_seq
-
-
 def _fetch_seq(loc: pd.DataFrame, src_fasta: Path) -> pd.Series:
-    """Fetch sequences from a source FASTA file, based on their record id, start and stop coordinates. 
+    """Fetch sequences from a source FASTA file, based on their record id, start and stop coordinates.
 
     Args:
-        loc (pd.DataFrame): A group of sequences in the same assembly. Required columns: 'record_idx', 'start' and 'stop'. 
-        src_fasta (Path): Path to the assembly FASTA file. 
+        loc (pd.DataFrame): A group of sequences in the same assembly. Required columns: 'record_idx', 'start' and 'stop'.
+        src_fasta (Path): Path to the assembly FASTA file.
 
     Returns:
-        pd.Series: Fetched sequences with the same index of `loc`. 
+        pd.Series: Fetched sequences with the same index of `loc`.
     """
     src_seq = load_fasta(src_fasta)
     # NOTE: assume all forward strand
     return loc.apply(
-        lambda row: src_seq[row['record_idx']][row['start']:row['stop']], 
+        lambda row: src_seq[row['record_idx']][row['start']:row['stop']],
         axis=1
     )
 
 
 def _get_paths_dl(taxa_list: list[str], prefix: Path, config: Config) -> list[Path]:
-    """Download assembly files for each taxon, and return the file paths. 
+    """Download assembly files for each taxon, and return the file paths.
 
     Args:
-        taxa_list (list[str]): See `tar_taxa` and `neg_taxa` in `Config` in `config.py`. 
-        prefix (Path): Download prefix. 
-        config (Config): See `Config` in `config.py`. 
+        taxa_list (list[str]): See `tar_taxa` and `neg_taxa` in `Config` in `config.py`.
+        prefix (Path): Download prefix.
+        config (Config): See `Config` in `config.py`.
 
     Returns:
-        list[Path]: Paths to assembly files. 
+        list[Path]: Paths to assembly files.
     """
     paths = list()
     # download genome assemblies under each taxon
     for taxon in taxa_list:
         download_paths = download_taxon(
-            taxon=taxon, 
-            prefix=prefix, 
-            level=config.level, 
-            source=config.source, 
-            annotated=config.annotated, 
-            exclude_mag=config.exclude_mag, 
-            gzip=config.gzip, 
+            taxon=taxon,
+            prefix=prefix,
+            level=config.level,
+            source=config.source,
+            annotated=config.annotated,
+            exclude_mag=config.exclude_mag,
+            gzip=config.gzip,
             api_key=config.api_key.get_secret_value() if config.api_key is not None else None,
-            overwrite=config.overwrite, 
+            overwrite=config.overwrite,
             n_cpu=config.n_cpu
         )
         if download_paths is not None:
@@ -431,13 +329,13 @@ def _get_paths_dl(taxa_list: list[str], prefix: Path, config: Config) -> list[Pa
 
 
 def _get_paths_txt(paths_txt: Path) -> list[Path]:
-    """Load assembly paths from a text file. 
+    """Load assembly paths from a text file.
 
     Args:
-        paths_txt (Path): See `tar_paths` and `neg_paths` in `Config` in `config.py`. 
+        paths_txt (Path): See `tar_paths` and `neg_paths` in `Config` in `config.py`.
 
     Returns:
-        list[Path]: Paths to assembly files. 
+        list[Path]: Paths to assembly files.
     """
     paths = load_paths_txt(paths_txt)
     logger.info(f'Found {len(paths)} assemblies from {paths_txt}')
@@ -445,13 +343,13 @@ def _get_paths_txt(paths_txt: Path) -> list[Path]:
 
 
 def _get_paths_dir(input_dir: Path) -> list[Path]:
-    """Load assembly paths from a directory (non-recursive). 
+    """Load assembly paths from a directory (non-recursive).
 
     Args:
-        input_dir (Path): See `tar_dir` and `neg_dir` in `Config` in `config.py`. 
+        input_dir (Path): See `tar_dir` and `neg_dir` in `Config` in `config.py`.
 
     Returns:
-        list[Path]: Paths to assembly files. 
+        list[Path]: Paths to assembly files.
     """
     paths = list()
 
@@ -470,16 +368,16 @@ def _get_paths_dir(input_dir: Path) -> list[Path]:
 
 
 def _download(config: Config, working_dir: Path) -> tuple[list[Path], list[Path]]:
-    """Download assemblies and return file paths. Return empty lists if nothing to download. 
+    """Download assemblies and return file paths. Return empty lists if nothing to download.
 
     Args:
-        config (Config): See `Config` in `config.py`. 
-        working_dir (Path): See `RunState` in `config.py`. 
+        config (Config): See `Config` in `config.py`.
+        working_dir (Path): See `RunState` in `config.py`.
 
     Returns:
         tuple: A tuple containing
-            1. list[Path]: Paths to downloaded target assemblies. 
-            2. list[Path]: Paths to downloaded non-target assemblies. 
+            1. list[Path]: Paths to downloaded target assemblies.
+            2. list[Path]: Paths to downloaded non-target assemblies.
     """
     tar_taxa = config.tar_taxa
     neg_taxa = config.neg_taxa
@@ -514,15 +412,15 @@ def _download(config: Config, working_dir: Path) -> tuple[list[Path], list[Path]
 
 
 def get_assemblies(config: Config, state: RunState) -> Assemblies:
-    """Load assembly paths and package them in an Assemblies instance. 
-    If taxonomy names are provided, download the genome FASTA files from NCBI. 
+    """Load assembly paths and package them in an Assemblies instance.
+    If taxonomy names are provided, download the genome FASTA files from NCBI.
 
     Args:
-        config (Config): See `Config` in `config.py`. 
-        state (RunState): See `RunState` in `config.py`. 
+        config (Config): See `Config` in `config.py`.
+        state (RunState): See `RunState` in `config.py`.
 
     Returns:
-        Assemblies: The Assemblies instance. 
+        Assemblies: The Assemblies instance.
     """
     tar_paths_txt = config.tar_paths
     neg_paths_txt = config.neg_paths
@@ -548,9 +446,9 @@ def get_assemblies(config: Config, state: RunState) -> Assemblies:
             neg_paths.extend(_get_paths_dir(neg_dir))
 
         if not tar_paths:
-            log_and_raise(RuntimeError, msg='No target assembly found.')
+            log_and_raise(RuntimeError, msg='No target assembly found')
         if not neg_paths:
-            log_and_raise(RuntimeError, msg='No non-target assembly found.')
+            log_and_raise(RuntimeError, msg='No non-target assembly found')
 
         # check if all paths are unique
         all_paths = tar_paths + neg_paths
diff --git a/src/seqwin/cli.py b/src/seqwin/cli.py
index b1bee51..4f7a19f 100644
--- a/src/seqwin/cli.py
+++ b/src/seqwin/cli.py
@@ -25,7 +25,7 @@
 
 
 app = typer.Typer(
-    help=f'Seqwin: Ultrafast identification of signature sequences', 
+    help=f'Seqwin: Ultrafast identification of signature sequences',
     add_completion=False, # do not show command completion options in help
     pretty_exceptions_show_locals=False, # do not show huge blocks of local variables
     add_help_option=False # disable built-in --help
@@ -48,158 +48,158 @@ def print_help(ctx: typer.Context, value: bool):
 def main(
     # even if only one value is provided, typer will still provide a list with only one element
     tar_taxa: list[str] | None = typer.Option(
-        None, '--tar-taxa', '-t', show_default=False, 
+        None, '--tar-taxa', '-t', show_default=False,
         help='Target NCBI taxonomy name or ID. Must be an exact match. '
-        'Repeat the option to pass multiple values (-t <tax1> -t <tax2> ...).', 
+        'Repeat the option to pass multiple values (-t <tax1> -t <tax2> ...).',
         rich_help_panel='Input selection'
-    ), 
+    ),
     neg_taxa: list[str] | None = typer.Option(
-        None, '--neg-taxa', '-n', show_default=False, 
+        None, '--neg-taxa', '-n', show_default=False,
         help='Non-target NCBI taxonomy name or ID. Must be an exact match. '
-        'Repeat the option to pass multiple values (-n <tax1> -n <tax2> ...).', 
+        'Repeat the option to pass multiple values (-n <tax1> -n <tax2> ...).',
         rich_help_panel='Input selection'
-    ), 
+    ),
     tar_paths: Path | None = typer.Option(
-        None, '--tar-paths', show_default=False, 
-        help='Text file containing paths to target genome FASTA files, one path per line. Gzipped FASTA is supported.', 
+        None, '--tar-paths', show_default=False,
+        help='Text file containing paths to target genome FASTA files, one path per line. Gzipped FASTA is supported.',
         rich_help_panel='Input selection'
-    ), 
+    ),
     neg_paths: Path | None = typer.Option(
-        None, '--neg-paths', show_default=False, 
-        help='Text file containing paths to non-target genome FASTA files, one path per line.', 
+        None, '--neg-paths', show_default=False,
+        help='Text file containing paths to non-target genome FASTA files, one path per line.',
         rich_help_panel='Input selection'
-    ), 
+    ),
     tar_dir: Path | None = typer.Option(
         None, '--tar-dir', show_default=False,
-        help='Directory containing target genome FASTA files.', 
+        help='Directory containing target genome FASTA files.',
         rich_help_panel='Input selection'
-    ), 
+    ),
     neg_dir: Path | None = typer.Option(
         None, '--neg-dir', show_default=False,
-        help='Directory containing non-target genome FASTA files.', 
+        help='Directory containing non-target genome FASTA files.',
         rich_help_panel='Input selection'
-    ), 
+    ),
     prefix: Path = typer.Option(
-        Path.cwd(), '--prefix', 
-        help='Parent path where the output directory will be created. Use the current working directory by default.', 
+        Path.cwd(), '--prefix',
+        help='Parent path where the output directory will be created. Use the current working directory by default.',
         rich_help_panel='Output options'
-    ), 
+    ),
     title: str = typer.Option(
-        'seqwin-out', '--title', '-o', 
-        help='Name of the output directory created under --prefix.', 
+        'seqwin-out', '--title', '-o',
+        help='Name of the output directory created under --prefix.',
         rich_help_panel='Output options'
-    ), 
+    ),
     overwrite: bool = typer.Option(
-        False, '--overwrite', show_default=False, 
-        help='Overwrite existing output files.', 
+        False, '--overwrite', show_default=False,
+        help='Overwrite existing output files.',
         rich_help_panel='Output options'
-    ), 
+    ),
     kmerlen: int = typer.Option(
-        21, '--kmerlen', '-k', 
-        help='K-mer length.', 
+        21, '--kmerlen', '-k',
+        help='K-mer length.',
         rich_help_panel='Signature options'
-    ), 
+    ),
     windowsize: int = typer.Option(
-        200, '--windowsize', '-w', 
-        help='Window size for minimizer sketch.', 
+        200, '--windowsize', '-w',
+        help='Window size for minimizer sketch.',
         rich_help_panel='Signature options'
-    ), 
+    ),
     penalty_th: float | None = typer.Option(
-        None, '--penalty-th', show_default=False, 
-        help='Node penalty threshold, from 0 to 1. If not provided, Seqwin computes it automatically.', 
+        None, '--penalty-th', show_default=False,
+        help='Node penalty threshold, from 0 to 1. If not provided, Seqwin computes it automatically.',
         rich_help_panel='Signature options'
-    ), 
+    ),
     # always default flags to False (can be reversed later in the function body)
     no_mash: bool = typer.Option(
-        False, '--no-mash', show_default=False, 
+        False, '--no-mash', show_default=False,
         help='Do not run Mash to estimate node penalty threshold. Instead, use minimizer sketches. '
         'This is much faster but the estimation might be biased. '
-        'Only used when --penalty-th is not provided.', 
+        'Only used when --penalty-th is not provided.',
         rich_help_panel='Signature options'
-    ), 
+    ),
     stringency: int = typer.Option(
-        5, '--stringency', '-s', show_default=True, 
+        5, '--stringency', '-s', show_default=True,
         help='Controls the sensitivity and specificity of output signatures (0-10). '
         'Increasing this value generally yields fewer and shorter signatures, while improving their sensitivity and specificity. '
         'Internally, Seqwin uses this setting to adjust the estimated node penalty threshold. '
-        'Only used when `--penalty-th` is not provided.', 
+        'Only used when `--penalty-th` is not provided.',
         rich_help_panel='Signature options'
-    ), 
+    ),
     min_len: int = typer.Option(
-        200, '--min-len', 
-        help='Minimum length of output signatures.', 
+        200, '--min-len',
+        help='Minimum length of output signatures.',
         rich_help_panel='Signature options'
-    ), 
+    ),
     max_len: int | None = typer.Option(
-        None, '--max-len', show_default=False, 
-        help='Estimated maximum length of output signatures. If not provided, no explicit limit is applied.', 
+        None, '--max-len', show_default=False,
+        help='Estimated maximum length of output signatures. If not provided, no explicit limit is applied.',
         rich_help_panel='Signature options'
-    ), 
+    ),
     no_blast: bool = typer.Option(
-        False, '--no-blast', show_default=False, 
-        help='Do not evaluate signature sequences with BLAST.', 
+        False, '--no-blast', show_default=False,
+        help='Do not evaluate signature sequences with BLAST.',
         rich_help_panel='Signature options'
-    ), 
+    ),
     # blast_neg_only: bool = typer.Option(
-    #     False, '--fast-blast', is_flag=True, flag_value=True, show_default=False, 
+    #     False, '--fast-blast', is_flag=True, flag_value=True, show_default=False,
     #     help='Only evaluate (BLAST) against non-target assemblies.'
-    # ), 
+    # ),
     level: Level = typer.Option(
         Level.contig, '--level', metavar='TEXT', # hide choices
         help="Limit downloads to genomes at or above this assembly level. "
-        "Possible values follow this order: 'contig', 'scaffold', 'chromosome', 'complete'", 
+        "Possible values follow this order: 'contig', 'scaffold', 'chromosome', 'complete'",
         rich_help_panel='NCBI download options'
-    ), 
+    ),
     source: Source = typer.Option(
         Source.genbank, '--source', metavar='TEXT', # hide choices
-        help="Genome source to download from. Supported values: 'genbank', 'refseq'", 
+        help="Genome source to download from. Supported values: 'genbank', 'refseq'",
         rich_help_panel='NCBI download options'
-    ), 
+    ),
     annotated: bool = typer.Option(
-        False, '--annotated', show_default=False, 
-        help='Only include annotated genomes.', 
+        False, '--annotated', show_default=False,
+        help='Only include annotated genomes.',
         rich_help_panel='NCBI download options'
-    ), 
+    ),
     exclude_mag: bool = typer.Option(
-        False, '--exclude-mag', show_default=False, 
-        help='Exclude metagenome-assembled genomes (MAGs).', 
+        False, '--exclude-mag', show_default=False,
+        help='Exclude metagenome-assembled genomes (MAGs).',
         rich_help_panel='NCBI download options'
-    ), 
+    ),
     no_gzip: bool = typer.Option(
-        False, '--no-gzip', show_default=False, 
-        help='Do not download genomes as gzipped FASTA.', 
+        False, '--no-gzip', show_default=False,
+        help='Do not download genomes as gzipped FASTA.',
         rich_help_panel='NCBI download options'
-    ), 
+    ),
     api_key: str | None = typer.Option(
         None, '--api-key', show_default=False,
-        help='NCBI API key passed to the Datasets command-line tools.', 
+        help='NCBI API key passed to the Datasets command-line tools.',
         rich_help_panel='NCBI download options'
-    ), 
+    ),
     download_only: bool = typer.Option(
-        False, '--download-only', show_default=False, 
-        help='Only download genome sequences without running Seqwin.', 
+        False, '--download-only', show_default=False,
+        help='Only download genome sequences without running Seqwin.',
         rich_help_panel='NCBI download options'
-    ), 
+    ),
     seed: int = typer.Option(
-        42, '--seed', 
-        help='Random seed for reproducibility.', 
+        42, '--seed',
+        help='Random seed for reproducibility.',
         rich_help_panel='Miscellaneous'
-    ), 
+    ),
     n_cpu: int = typer.Option(
-        4, '--threads', '-p', 
-        help='Number of parallel processes or threads to use.', 
+        4, '--threads', '-p',
+        help='Number of parallel processes or threads to use.',
         rich_help_panel='Miscellaneous'
-    ), 
+    ),
     version: bool = typer.Option(
-        False, '--version', callback=print_version, show_default=False, expose_value=False, 
+        False, '--version', callback=print_version, show_default=False, expose_value=False,
         is_eager=True, # run this before any other options
-        help='Show Seqwin version and exit.', 
+        help='Show Seqwin version and exit.',
         rich_help_panel='Miscellaneous'
-    ), 
+    ),
     help_: bool = typer.Option(
-        False, '--help', '-h', callback=print_help, show_default=False, expose_value=False, 
+        False, '--help', '-h', callback=print_help, show_default=False, expose_value=False,
         is_eager=True, # run this before any other options
-        help='Show this message and exit.', 
+        help='Show this message and exit.',
         rich_help_panel='Miscellaneous'
     )
 ):
@@ -210,31 +210,31 @@ def main(
             raise typer.BadParameter('You must provide at least one non-target input: --neg-paths, --neg-taxa, or --neg-dir')
 
     config = Config(
-        tar_taxa=tar_taxa, 
-        neg_taxa=neg_taxa, 
-        tar_paths=tar_paths, 
-        neg_paths=neg_paths, 
+        tar_taxa=tar_taxa,
+        neg_taxa=neg_taxa,
+        tar_paths=tar_paths,
+        neg_paths=neg_paths,
         tar_dir=tar_dir,
         neg_dir=neg_dir,
-        prefix=prefix, 
-        title=title, 
-        overwrite=overwrite, 
-        kmerlen=kmerlen, 
-        windowsize=windowsize, 
-        penalty_th=penalty_th, 
-        run_mash=not no_mash, 
-        stringency=stringency, 
-        min_len=min_len, 
-        max_len=max_len, 
-        run_blast=not no_blast, 
-        #blast_neg_only=blast_neg_only, 
-        seed=seed, 
-        n_cpu=n_cpu, 
-        level=level, 
-        source=source, 
-        annotated=annotated, 
-        exclude_mag=exclude_mag, 
-        gzip=not no_gzip, 
+        prefix=prefix,
+        title=title,
+        overwrite=overwrite,
+        kmerlen=kmerlen,
+        windowsize=windowsize,
+        penalty_th=penalty_th,
+        run_mash=not no_mash,
+        stringency=stringency,
+        min_len=min_len,
+        max_len=max_len,
+        run_blast=not no_blast,
+        #blast_neg_only=blast_neg_only,
+        seed=seed,
+        n_cpu=n_cpu,
+        level=level,
+        source=source,
+        annotated=annotated,
+        exclude_mag=exclude_mag,
+        gzip=not no_gzip,
         api_key=api_key,
         download_only=download_only
     )
diff --git a/src/seqwin/config.py b/src/seqwin/config.py
index a99898b..1ea5874 100644
--- a/src/seqwin/config.py
+++ b/src/seqwin/config.py
@@ -2,7 +2,7 @@
 Configurations
 ==============
 
-Seqwin run configurations. Including user/dev configs and internal configs. 
+Seqwin run configurations. Including user/dev configs and internal configs.
 
 Dependencies:
 -------------
@@ -45,10 +45,10 @@
 
 # init root logger
 logging.basicConfig(
-    format=_LOG_FMT, 
-    datefmt=_LOG_DATEFMT, 
-    level=logging.INFO, 
-    stream=sys.stdout, 
+    format=_LOG_FMT,
+    datefmt=_LOG_DATEFMT,
+    level=logging.INFO,
+    stream=sys.stdout,
 )
 
 from pathlib import Path
@@ -72,7 +72,7 @@
 
 
 class Config(BaseModel):
-    """Seqwin configurations. 
+    """Seqwin configurations.
 
     Attributes:
         tar_taxa (list[str] | None): Target NCBI taxonomy name(s) / ID(s). Must be exact matches. [None]
@@ -114,7 +114,7 @@ class Config(BaseModel):
 
         seed (int): Random seed for reproducibility. [42]
         n_cpu (int): Number of parallel processes or threads to use. [4]
-        version (str): Seqwin version. 
+        version (str): Seqwin version.
     """
     # Inputs
     tar_taxa: list[str] | None = None
@@ -192,7 +192,7 @@ def _check_inputs(self) -> 'Config':
         if (not HAS_DATASETS) and (self.tar_taxa or self.neg_taxa):
             raise FileNotFoundError(
                 ('ncbi-datasets-cli is not installed. Genomes cannot be downloaded from the '
-                'provided taxon names or IDs. Please provide local file paths instead.'))
+                'provided taxon names or IDs. Please provide local files instead'))
 
         if not self.download_only:
             if (self.tar_paths is None) and (self.tar_taxa is None) and (self.tar_dir is None):
@@ -213,27 +213,27 @@ def _check_inputs(self) -> 'Config':
 
     model_config = {
         # similar to @dataclass(slots=True, frozen=True)
-        'frozen': True, 
-        'slots': True, 
-        'validate_default': True, 
+        'frozen': True,
+        'slots': True,
+        'validate_default': True,
         'hide_input_in_errors': True
     }
 
 
 @dataclass(slots=True)
 class RunState:
-    """Runtime variables. 
+    """Runtime variables.
 
     Attributes:
-        working_dir (Path): Working directory, defined by prefix and title. 
-        rng (random.Random): Built-in `random.Random` instance for reproducibility. 
-        n_tar (int | None): Number of target assemblies. 
-        n_neg (int | None): Number of non-target assemblies. 
-        penalty_th (float | None): Node penalty threshold (user input or auto-computed). 
-        edge_weight_th (float | None): Graph edge weight threshold. 
-        min_nodes (int | None): Min number of nodes for a low-penalty subgraph. 
-        max_nodes (int | None): Max number of nodes for a low-penalty subgraph. 
-        blastdb (Path | None): Path to the BLAST database inside the working directory. 
+        working_dir (Path): Working directory, defined by prefix and title.
+        rng (random.Random): Built-in `random.Random` instance for reproducibility.
+        n_tar (int | None): Number of target assemblies.
+        n_neg (int | None): Number of non-target assemblies.
+        penalty_th (float | None): Node penalty threshold (user input or auto-computed).
+        edge_weight_th (float | None): Graph edge weight threshold.
+        min_nodes (int | None): Min number of nodes for a low-penalty subgraph.
+        max_nodes (int | None): Max number of nodes for a low-penalty subgraph.
+        blastdb (Path | None): Path to the BLAST database inside the working directory.
     """
     working_dir: Path
     rng: Random
@@ -248,7 +248,7 @@ class RunState:
 
 @dataclass(slots=True, frozen=True)
 class WorkingDir:
-    """Files and directories under the working directory. 
+    """Files and directories under the working directory.
 
     Attributes:
         log (str): Seqwin log file. ['seqwin.log']
@@ -276,7 +276,7 @@ class WorkingDir:
 
 @dataclass(slots=True, frozen=True)
 class BlastConfig:
-    """Settings for the BLAST commands `makeblastdb` and `blastn`. 
+    """Settings for the BLAST commands `makeblastdb` and `blastn`.
 
     Attributes:
         title_neg_only (str): DB title when created from non-target assemblies only. ['neg-only']
@@ -286,7 +286,7 @@ class BlastConfig:
         str2bool (Mapping[str, bool]): Reversed mapping of bool2str for parsing BLAST outputs. ['y' -> True, 'n' -> False]
         header_sep (str): Separator used in FASTA headers. Should pick a rare char (cannot be '$', BLAST treats it as a special char). ['@']
         task (str): Presets for BLAST parameters ('blastn', 'blastn-short', 'megablast'). ['blastn']
-        columns (tuple[str, ...]): Columns to be included in the BLAST TSV output. See `ncbi.py` for more information. 
+        columns (tuple[str, ...]): Columns to be included in the BLAST TSV output. See `ncbi.py` for more information.
         batch_size (int): Number of query sequences in a single BLAST run. [1000]
     """
     title_neg_only: str = 'neg-only'
@@ -302,33 +302,33 @@ class BlastConfig:
     header_sep: str = '@'
     task: Task = Task.blastn
     columns = (
-        'qseqid', 
-        'sseqid', 
-        'nident', 
-        'mismatch', 
-        'gaps', 
-        'qstart', 
-        'qend', 
-        'sstart', 
-        'send', 
-        'evalue', 
-        'bitscore', 
+        'qseqid',
+        'sseqid',
+        'nident',
+        'mismatch',
+        'gaps',
+        'qstart',
+        'qend',
+        'sstart',
+        'send',
+        'evalue',
+        'bitscore',
         'sseq'
     )
     batch_size: int = 1000
 
 
 def config_logger(file: Path, level: int) -> None:
-    """Add a file handler and set logging level for the root logger. 
+    """Add a file handler and set logging level for the root logger.
 
     Args:
-        file (Path): Path to the log file. 
-        level (int): Logging level (e.g., `logging.INFO`). 
+        file (Path): Path to the log file.
+        level (int): Logging level (e.g., `logging.INFO`).
     """
     # use the same format
     logFormatter = logging.Formatter(
-        fmt=_LOG_FMT, 
-        datefmt=_LOG_DATEFMT, 
+        fmt=_LOG_FMT,
+        datefmt=_LOG_DATEFMT,
         style='%'
     )
     fileHandler = logging.FileHandler(file, mode='a')
diff --git a/src/seqwin/core.py b/src/seqwin/core.py
index d897c18..9a7d8be 100644
--- a/src/seqwin/core.py
+++ b/src/seqwin/core.py
@@ -2,7 +2,7 @@
 Core
 ====
 
-Seqwin entry point. 
+Seqwin entry point.
 
 Dependencies:
 -------------
@@ -42,15 +42,15 @@
 
 
 class Seqwin(object):
-    """Seqwin run instance. 
+    """Seqwin run instance.
 
     Attributes:
-        config (Config): See `Config` in `config.py`. 
-        state (RunState): See `RunState` in `config.py`. 
-        assemblies (Assemblies): See `Assemblies` in `assemblies.py`. 
-        kmers (KmerGraph | None): See `KmerGraph` in `kmers.py`. Generated with `self.run()`. 
-        mash (pd.DataFrame | None): Tabular output of `mash dist`. Generated with `self.run()`. 
-        markers (list[ConnectedKmers] | None): See `ConnectedKmers` in `markers.py`. Generated with `self.run()`. 
+        config (Config): See `Config` in `config.py`.
+        state (RunState): See `RunState` in `config.py`.
+        assemblies (Assemblies): See `Assemblies` in `assemblies.py`.
+        kmers (KmerGraph | None): See `KmerGraph` in `kmers.py`. Generated with `self.run()`.
+        mash (pd.DataFrame | None): Tabular output of `mash dist`. Generated with `self.run()`.
+        markers (list[ConnectedKmers] | None): See `ConnectedKmers` in `markers.py`. Generated with `self.run()`.
     """
     __slots__ = ('config', 'state', 'assemblies', 'kmers', 'mash', 'markers')
     config: Config
@@ -61,14 +61,14 @@ class Seqwin(object):
     markers: list[ConnectedKmers] | None
 
     def __init__(self, config: Config) -> None:
-        """Initiate a Seqwin run instance. 
-        1. Create a working directory. 
-        2. Initialize the logger. 
-        3. Save config to JSON. 
-        4. Load all assemblies. 
+        """Initiate a Seqwin run instance.
+        1. Create a working directory.
+        2. Initialize the logger.
+        3. Save config to JSON.
+        4. Load all assemblies.
 
         Args:
-            config (Config): See `Config` in `config.py`. 
+            config (Config): See `Config` in `config.py`.
         """
         prefix = config.prefix
         title = config.title
@@ -119,7 +119,7 @@ def __init__(self, config: Config) -> None:
         self.markers = None
 
     def run(self) -> None:
-        """Build the k-mer graph and extract candidate markers. 
+        """Build the k-mer graph and extract candidate markers.
         """
         config = self.config
         state = self.state
@@ -148,13 +148,13 @@ def run(self) -> None:
 
 
 def run(config: Config) -> Seqwin:
-    """Run Seqwin. 
+    """Run Seqwin.
 
     Args:
-        config (Config): See `Config` in `config.py`. 
+        config (Config): See `Config` in `config.py`.
 
     Returns:
-        Seqwin: The Seqwin run instance. 
+        Seqwin: The Seqwin run instance.
     """
     seqwin = Seqwin(config)
     if not config.download_only:
@@ -163,13 +163,13 @@ def run(config: Config) -> Seqwin:
 
 
 def load(path: str | Path) -> Seqwin:
-    """Load a Seqwin run instance from file. 
+    """Load a Seqwin run instance from file.
 
     Args:
-        path (str | Path): Path to the Seqwin run snapshot (`results.seqwin`). 
+        path (str | Path): Path to the Seqwin run snapshot (`results.seqwin`).
 
     Returns:
-        Seqwin: The Seqwin run instance. 
+        Seqwin: The Seqwin run instance.
     """
     if isinstance(path, str):
         path = Path(path)
diff --git a/src/seqwin/graph/__init__.py b/src/seqwin/graph/__init__.py
index 721c1d1..fdca009 100644
--- a/src/seqwin/graph/__init__.py
+++ b/src/seqwin/graph/__init__.py
@@ -2,7 +2,7 @@
 Minimizer Graph
 ===============
 
-Core functions and dtypes for Seqwin minimizer graphs. 
+Core functions and dtypes for Seqwin minimizer graphs.
 
 Usage:
 ------
@@ -10,11 +10,11 @@
 >>> from pathlib import Path
 >>> from seqwin.graph import build
 >>> kmers, idx, nodes, edges, record_ids = build(
->>>     assembly_paths=[Path('example1.fa'), Path('example2.fa.gz')], 
->>>     kmerlen=21, 
->>>     windowsize=200, 
->>>     assembly_idx=[0, 1], 
->>>     is_target=[True, False], 
+>>>     assembly_paths=[Path('example1.fa'), Path('example2.fa.gz')],
+>>>     kmerlen=21,
+>>>     windowsize=200,
+>>>     assembly_idx=[0, 1],
+>>>     is_target=[True, False],
 >>> )
 ```
 
@@ -45,17 +45,17 @@
 from .utils import OrderedKmers
 
 KMER_DTYPE = np.dtype([
-    ('pos', np.uint32), 
-    ('record_idx', np.uint16), 
+    ('pos', np.uint32),
+    ('record_idx', np.uint16),
     ('assembly_idx', np.uint16)
 ])
 
 NODE_DTYPE = np.dtype([
-    ('hash', np.uint64), 
-    ('start', np.uint64), 
-    ('stop', np.uint64), 
-    ('n_tar', np.uint32), 
-    ('n_neg', np.uint32), 
+    ('hash', np.uint64),
+    ('start', np.uint64),
+    ('stop', np.uint64),
+    ('n_tar', np.uint32),
+    ('n_neg', np.uint32),
     ('penalty', np.float64)
 ])
 
@@ -67,23 +67,39 @@
 
 
 def build(
-    assembly_paths: Iterable[Path], 
-    kmerlen: int, 
-    windowsize: int, 
-    assembly_idx: Iterable[int], 
-    is_target: Iterable[bool], 
+    assembly_paths: Iterable[Path],
+    kmerlen: int,
+    windowsize: int,
+    assembly_indices: Iterable[int],
+    is_targets: Iterable[bool],
     n_cpu: int = 1
 ) -> tuple[
-    NDArray[np.void], 
-    NDArray[np.uint64], 
-    NDArray[np.void], 
-    NDArray[np.void], 
+    NDArray[np.void],
+    NDArray[np.uint64],
+    NDArray[np.void],
+    NDArray[np.void],
     list[tuple[str, ...]]
 ]:
-    """Build a Seqwin minimizer graph. 
-    - `assembly_paths`, `assembly_idx`, and `is_target` are parallel lists. 
-    - In the returned arrays, `kmers` is grouped by node/hash. 
-    - For every node:
+    """Build a Seqwin minimizer graph.
+
+    Example usage:
+    ```python
+    >>> from seqwin.graph import build
+    >>> kmers, idx, nodes, edges, record_ids = build(
+    >>>     assembly_paths = ...,
+    >>>     kmerlen = 21,
+    >>>     windowsize = 200,
+    >>>     assembly_indices = ...,
+    >>>     is_targets = ...,
+    >>>     n_cpu = 4,
+    >>> )
+    ```
+    - `assembly_paths`, `assembly_indices`, and `is_targets` are parallel lists.
+    - `kmers` stores minimizer occurrences in all assemblies, grouped and sorted by hash.
+    - `idx` is parallel to `kmers` and stores each minimizer's original generation index, ordered by genomic position.
+    - `nodes` and `edges` are sorted by hash.
+
+    The `[start, stop)` range in each node identifies minimizers with this hash.
     ```python
     >>> kmer_group = kmers[node['start']:node['stop']]
     >>> group_hash = node['hash']
@@ -91,60 +107,70 @@ def build(
     ```
 
     Args:
-        assembly_paths (Iterable[Path]): Path to each assembly file in FASTA format (gzip supported). 
-        kmerlen (int): k-mer length. 
-        windowsize (int): Window size for minimizer sketch. 
-        assembly_idx (Iterable[int]): Index of each assembly. 
-        is_target (Iterable[bool]): True for target assemblies. 
-        n_cpu (int, optional): Number of threads. [1]
+        assembly_paths (Iterable[Path]): Paths to input assemblies in FASTA format (plain or gzipped).
+        kmerlen (int): K-mer length for minimizer sketch.
+        windowsize (int): Window size for minimizer sketch.
+        assembly_indices (Iterable[int]): Assembly indices.
+        is_targets (Iterable[bool]): Whether each assembly is a target assembly.
+        n_cpu (int, optional): Number of worker threads to use. [1]
 
     Returns:
         tuple: A tuple containing
-            1. NDArray[np.void]: A 1-D NumPy structured array of k-mers from all assemblies, with dtype `KMER_DTYPE`. 
-                Each element represents a minimizer, with fields, 
-                - 'pos' (uint32): Position of the first base of the minimizer. 
-                - 'record_idx' (uint16): 0-based index of the sequence records, in the same order as they appear in the FASTA file. 
-                - 'assembly_idx' (uint16): Assembly index. 
-            2. NDArray[np.uint64]: The original indices assigned when k-mers are generated (ordered by genomic positions). 
-            3. NDArray[np.void]: A 1-D NumPy structured array of k-mer nodes, with dtype `NODE_DTYPE`. 
-            4. NDArray[np.void]: A 1-D NumPy structured array of weighted, undirected edges, with dtype `EDGE_DTYPE`. 
-            5. list[tuple[str, ...]]: FASTA record IDs of each assembly. 
+            1. NDArray[np.void]: A 1-D NumPy structured array of minimizers from all assemblies.
+                Dtype: `KMER_DTYPE`
+                - 'pos' (uint32): 0-based position of the minimizer within its FASTA record.
+                - 'record_idx' (uint16): 0-based index of the FASTA record within the assembly.
+                - 'assembly_idx' (uint16): Assembly index assigned to the source assembly.
+            2. NDArray[np.uint64]: The original indices assigned when minimizers are generated (ordered by genomic positions).
+            3. NDArray[np.void]: A 1-D NumPy structured array of minimizer nodes.
+                Dtype: `NODE_DTYPE`
+                - 'hash' (uint64): Hash value of the minimizers represented by this node.
+                - 'start' (uint64): Start of the half-open range for this node's minimizer entries.
+                - 'stop' (uint64): End of the half-open range for this node's minimizer entries.
+                - 'n_tar' (uint32): Number of target assemblies containing this minimizer hash.
+                - 'n_neg' (uint32): Number of non-target assemblies containing this minimizer hash.
+                - 'penalty' (float64): Node penalty score used for downstream graph filtering.
+            4. NDArray[np.void]: A 1-D NumPy structured array of weighted, undirected edges.
+                Dtype: `EDGE_DTYPE`
+                - 'first' (uint64): Smaller endpoint hash of the undirected edge.
+                - 'second' (uint64): Larger endpoint hash of the undirected edge.
+                - 'weight' (uint64): Number of assemblies where the endpoints are adjacent.
+            5. list[tuple[str, ...]]: FASTA record IDs of each assembly.
     """
-    kmers, idx, nodes, edges, idx_to_id = _build_native(
-        list(str(p) for p in assembly_paths), 
-        int(kmerlen), 
-        int(windowsize), 
-        list(int(idx) for idx in assembly_idx), 
-        list(bool(target) for target in is_target), 
+    return _build_native(
+        list(str(p) for p in assembly_paths),
+        int(kmerlen),
+        int(windowsize),
+        list(int(i) for i in assembly_indices),
+        list(bool(t) for t in is_targets),
         int(n_cpu)
     )
-    return kmers, idx, nodes, edges, idx_to_id
 
 
 def _filter_kmers(
-    kmers: NDArray[np.void], 
-    idx: NDArray[np.uint64], 
-    nodes: NDArray[np.void], 
+    kmers: NDArray[np.void],
+    idx: NDArray[np.uint64],
+    nodes: NDArray[np.void],
     used_hashes: frozenset[np.uint64]
 ) -> tuple[
-    NDArray[np.void], 
-    NDArray[np.uint64], 
+    NDArray[np.void],
+    NDArray[np.uint64],
     NDArray[np.void]
 ]:
     """
-    1. Remove k-mers (`kmers`, `idx` and `nodes`) not included in `used_hashes`. 
-    2. Update 'start' and 'stop' in nodes. 
+    1. Remove k-mers (`kmers`, `idx` and `nodes`) not included in `used_hashes`.
+    2. Update 'start' and 'stop' in nodes.
 
     Args:
-        kmers (NDArray): See `KmerGraph.kmers` (already sorted by 'hash'). 
-        idx (NDArray): See `KmerGraph.idx`. 
-        nodes (NDArray): See `KmerGraph.nodes` (sorted by 'hash'). 
-        used_hashes (frozenset[np.uint64]): K-mers and nodes with these hash values are kept. 
+        kmers (NDArray): See `KmerGraph.kmers`.
+        idx (NDArray): See `KmerGraph.idx`.
+        nodes (NDArray): See `KmerGraph.nodes`.
+        used_hashes (frozenset[np.uint64]): K-mers and nodes with these hash values are kept.
 
     Returns:
         tuple: A tuple containing
-            1. NDArray: See `KmerGraph.kmers`. 
-            2. NDArray: See `KmerGraph.idx`. 
-            3. NDArray: See `KmerGraph.nodes`. 
+            1. NDArray: See `KmerGraph.kmers`.
+            2. NDArray: See `KmerGraph.idx`.
+            3. NDArray: See `KmerGraph.nodes`.
     """
     return _filter_kmers_native(kmers, idx, nodes, used_hashes)
diff --git a/src/seqwin/graph/utils.py b/src/seqwin/graph/utils.py
index 681c442..b6d0856 100644
--- a/src/seqwin/graph/utils.py
+++ b/src/seqwin/graph/utils.py
@@ -2,7 +2,7 @@
 Graph
 =====
 
-Graph utilities. 
+Graph utilities.
 
 Dependencies:
 -------------
@@ -47,25 +47,25 @@
 
 
 class WeightedGraph(Counter):
-    """A weighted graph object inherited from `collections.Counter`, with data structure of {edge: weight}. 
-    An edge should contain two hashable elements and ordering matters, with edge direction of (first -> second). 
+    """A weighted graph object inherited from `collections.Counter`, with data structure of {edge: weight}.
+    An edge should contain two hashable elements and ordering matters, with edge direction of (first -> second).
     """
     def __init__(self, edges: Iterable[tuple]=()) -> None:
-        """Create a weighted graph with an Iterable of edges. 
-        Each edge should be a tuple containing two hashable elements as nodes. 
-        Note that the ordering of the two nodes matters, with edge direction of (first node -> second node). 
+        """Create a weighted graph with an Iterable of edges.
+        Each edge should be a tuple containing two hashable elements as nodes.
+        Note that the ordering of the two nodes matters, with edge direction of (first node -> second node).
 
         Args:
-            edges (Iterable[tuple], optional): Each edge should be a tuple containing two hashable elements as nodes. 
-            Note that the ordering of the two nodes matters, with edge direction of (first node -> second node). 
+            edges (Iterable[tuple], optional): Each edge should be a tuple containing two hashable elements as nodes.
+            Note that the ordering of the two nodes matters, with edge direction of (first node -> second node).
         """
         super().__init__(edges)
 
     def add_path(self, nodes: Iterable, cyclic=False) -> None:
-        """Add a path to the weighted graph. 
+        """Add a path to the weighted graph.
 
         Args:
-            nodes (Iterable): A path will be constructed from the nodes (in order) and added to the graph. 
+            nodes (Iterable): A path will be constructed from the nodes (in order) and added to the graph.
             cyclic (bool, optional): If True, connect the last node back to the first node. [False]
         """
         nodes = iter(nodes)
@@ -77,38 +77,38 @@ def add_path(self, nodes: Iterable, cyclic=False) -> None:
             return
         if cyclic:
             stop_nodes = chain(stop_nodes, (first_node,))
-        
+
         # add edges to graph
         self.update(
             tuple((u, v)) for u, v in zip(start_nodes, stop_nodes)
         )
 
     def to_nxGraph(self) -> nx.Graph:
-        """Convert to networkx.Graph (undirected), with edge weights set as edge attribute EDGE_W. 
+        """Convert to networkx.Graph (undirected), with edge weights set as edge attribute EDGE_W.
         """
         return nx.Graph((*edge, {EDGE_W: weight}) for edge, weight in self.items())
 
 
 class OrderedKmers(tuple):
-    """Ordered k-mers created from an Iterable of integers. 
-    The `which_strand()` method can take another Iterable of k-mers and determine its strand ('+'/'-'/'?'/'u'), 
-    by comparing its ordering to self. 
+    """Ordered k-mers created from an Iterable of integers.
+    The `which_strand()` method can take another Iterable of k-mers and determine its strand ('+'/'-'/'?'/'u'),
+    by comparing its ordering to self.
 
     Attributes:
-        rev (tuple): K-mers in reversed order. 
-        is_dup (bool): True if there are duplicates in the k-mers. 
-        warning (set): For debugging only. 
+        rev (tuple): K-mers in reversed order.
+        is_dup (bool): True if there are duplicates in the k-mers.
+        warning (set): For debugging only.
 
     Examples:
         ```
         l = [
             (1,2,3,3,4,5),
             (5,4,3,3,2,1),
-            (1,2,3,4,5), 
-            (5,4,3,2,1), 
-            (2,), 
-            (0,), 
-            (6,5), 
+            (1,2,3,4,5),
+            (5,4,3,2,1),
+            (2,),
+            (0,),
+            (6,5),
             (9,10),
             (1,3,5),
             (2,3,4),
@@ -131,12 +131,12 @@ def __new__(cls, kmers: Iterable[int]):
         return super().__new__(cls, kmers)
 
     def __init__(self, kmers: Iterable[int]) -> None:
-        """Ordered k-mers created from an Iterable of integers. 
-        The `which_strand()` method can take another Iterable of k-mers and determine its strand ('+'/'-'/'?'/'u'), 
-        by comparing its ordering to self. 
+        """Ordered k-mers created from an Iterable of integers.
+        The `which_strand()` method can take another Iterable of k-mers and determine its strand ('+'/'-'/'?'/'u'),
+        by comparing its ordering to self.
 
         Args:
-            kmers (Iterable[int]): K-mers as an Iterable of integers. 
+            kmers (Iterable[int]): K-mers as an Iterable of integers.
         """
         # here self is already created as a tuple
         # kmers is not used here, but have to keep it or it will raise a TypeError (for docstring as well)
@@ -144,19 +144,19 @@ def __init__(self, kmers: Iterable[int]) -> None:
         self._idx_map = {kmer: idx for idx, kmer in enumerate(self)}
         self.is_dup = len(self._idx_map) < self.__len__() # True if there are duplicated k-mers
         self.warning = set()
-    
+
     def which_strand(self, kmers: Iterable[int]) -> str:
-        """Given an Iterable of k-mers, compare its ordering to self and determine its strand ('+'/'-'/'?'/'u'). 
+        """Given an Iterable of k-mers, compare its ordering to self and determine its strand ('+'/'-'/'?'/'u').
 
         Args:
-            kmers (Iterable[int]): K-mers as an Iterable of integers. 
+            kmers (Iterable[int]): K-mers as an Iterable of integers.
 
         Returns:
             str: strand type
-            - '+': forward strand, 
-            - '-': reverse strand, 
-            - '?': unknown strand, 
-            - 'u': only one shared k-mer with self, so the strand has to be determined by other methods. 
+            - '+': forward strand,
+            - '-': reverse strand,
+            - '?': unknown strand,
+            - 'u': only one shared k-mer with self, so the strand has to be determined by other methods.
         """
         # keep in mind that there might be k-mers not found in self
         idx_map = self._idx_map
@@ -223,20 +223,20 @@ def check_order(orderedKmers) -> bool:
 
 
 def compose_weighted_graphs(graphs: Iterable[WeightedGraph]):
-    """Compose multiple WeightedGraph objects by adding the weights of same edges together. 
+    """Compose multiple WeightedGraph objects by adding the weights of same edges together.
 
     Args:
-        graphs (Iterable[WeightedGraph]): Graphs to be composed. 
+        graphs (Iterable[WeightedGraph]): Graphs to be composed.
 
     Returns:
-        WeightedGraph: The composed graph. 
+        WeightedGraph: The composed graph.
     """
     graphs = iter(graphs)
     try:
         merged_graph = next(graphs)
     except StopIteration:
         raise ValueError('No graph is given to compose. ')
-    
+
     merged_graph = merged_graph.copy() # a shallow copy
     for g in graphs:
         merged_graph.update(g)
@@ -244,11 +244,11 @@ def compose_weighted_graphs(graphs: Iterable[WeightedGraph]):
 
 
 def add_path_weighted(graph: nx.Graph, path: Sequence) -> None:
-    """Add a path to a weighted, undirected graph. Increment edge weight by 1 if the edge already exists. 
+    """Add a path to a weighted, undirected graph. Increment edge weight by 1 if the edge already exists.
 
     Args:
-        graph (nx.Graph): A weighted, undirected graph. 
-        path (Sequence): A path (sequence of nodes) to be added to graph. 
+        graph (nx.Graph): A weighted, undirected graph.
+        path (Sequence): A path (sequence of nodes) to be added to graph.
     """
     # loop through each consecutive pair of nodes in the path
     for i in range(len(path) - 1):
@@ -264,19 +264,19 @@ def add_path_weighted(graph: nx.Graph, path: Sequence) -> None:
 
 if _HAS_MPL:
     def draw_weighted_graph(
-        graph: nx.Graph, 
-        save_path: str | None=None, 
-        figsize: tuple | None=None, 
-        node_size: int=200, 
-        edge_width: int=2, 
-        font_size: int=8, 
+        graph: nx.Graph,
+        save_path: str | None=None,
+        figsize: tuple | None=None,
+        node_size: int=200,
+        edge_width: int=2,
+        font_size: int=8,
         seed: int=0
     ) -> None:
-        """Draw a NetworkX graph with edge attribute 'weight'. 
+        """Draw a NetworkX graph with edge attribute 'weight'.
         Code adapted from `networkx doc<https://networkx.org/documentation/stable/auto_examples/drawing/plot_weighted_graph.html>`__.
 
         Args:
-            graph (nx.Graph): A weighted, undirected graph. 
+            graph (nx.Graph): A weighted, undirected graph.
             save_path (str | None, optional): Path to save the figure in SVG format. None for showing the figure without saving. [None]
         """
         # positions for all nodes - seed for reproducibility
@@ -284,7 +284,7 @@ def draw_weighted_graph(
 
         if figsize is not None:
             plt.figure(figsize=figsize)
-        
+
         # nodes
         nx.draw_networkx_nodes(graph, pos, node_size=node_size)
 
diff --git a/src/seqwin/helpers.py b/src/seqwin/helpers.py
index b220a5f..ebaf844 100644
--- a/src/seqwin/helpers.py
+++ b/src/seqwin/helpers.py
@@ -2,7 +2,7 @@
 Helpers
 =======
 
-Helper functions for `kmers.py`. 
+Helper functions for `kmers.py`.
 
 Dependencies:
 -------------
@@ -33,39 +33,39 @@
 
 
 def get_subgraphs(
-    graph: nx.Graph, 
-    penalty_th: float, 
-    min_nodes: int, 
-    max_nodes: int | None, 
+    graph: nx.Graph,
+    penalty_th: float,
+    min_nodes: int,
+    max_nodes: int | None,
     rng: Random
 ) -> tuple[
-    tuple[frozenset[np.uint64], ...], 
+    tuple[frozenset[np.uint64], ...],
     frozenset[np.uint64]
 ]:
-    """Find disjoint (no shared node) subgraphs whose average node-penalty ≤ `penalty_th` and size within `size_th`. 
-    1. Remove low-weight edges and isolated nodes from `graph`. 
-    2. Find nodes with penalty ≤ `penalty_th` as seeds of subgraphs. 
-    3. Greedy seed-expansion with breadth first search (BFS), where the neighboring node with the lowest penalty is 
-        selected in each iteration. 
-
-    A heap frontier (nodes to be visited in BFS) is used to accelerate the expansion process. 
-    The heap is implemented with the built-in Python `heapq` module, which is a min-heap. 
-    E.g., when tuples of `(penalty, node)` are pushed to the heap, it will always pop the tuple with the smallest `penalty` first. 
-    This is faster than calling `min()` every time to fetch the node with the lowest penalty. 
-    When tested on the Salmonella dataset (576 genomes, no edge filtering), this implementation is more than 3x faster than the naive one. 
-    However, the performance gain becomes less significant when more low-weight edges are removed. 
+    """Find disjoint (no shared node) subgraphs whose average node-penalty ≤ `penalty_th` and size within `size_th`.
+    1. Remove low-weight edges and isolated nodes from `graph`.
+    2. Find nodes with penalty ≤ `penalty_th` as seeds of subgraphs.
+    3. Greedy seed-expansion with breadth first search (BFS), where the neighboring node with the lowest penalty is
+        selected in each iteration.
+
+    A heap frontier (nodes to be visited in BFS) is used to accelerate the expansion process.
+    The heap is implemented with the built-in Python `heapq` module, which is a min-heap.
+    E.g., when tuples of `(penalty, node)` are pushed to the heap, it will always pop the tuple with the smallest `penalty` first.
+    This is faster than calling `min()` every time to fetch the node with the lowest penalty.
+    When tested on the Salmonella dataset (576 genomes, no edge filtering), this implementation is more than 3x faster than the naive one.
+    However, the performance gain becomes less significant when more low-weight edges are removed.
 
     Args:
-        graph (nx.Graph): See `KmerGraph.graph`. 
-        penalty_th (float): See `Config` in `config.py`. 
-        min_nodes (int): See `Config` in `config.py`. 
-        max_nodes (int | None): See `Config` in `config.py`. 
-        rng (random.Random): See `RunState` in `config.py`. 
+        graph (nx.Graph): See `KmerGraph.graph`.
+        penalty_th (float): See `Config` in `config.py`.
+        min_nodes (int): See `Config` in `config.py`.
+        max_nodes (int | None): See `Config` in `config.py`.
+        rng (random.Random): See `RunState` in `config.py`.
 
     Returns:
         tuple: A tuple containing
-            1. tuple[frozenset[np.uint64], ...]: See `KmerGraph.subgraphs`. 
-            2. frozenset[np.uint64]: Union of k-mer hash values in all subgraphs. 
+            1. tuple[frozenset[np.uint64], ...]: See `KmerGraph.subgraphs`.
+            2. frozenset[np.uint64]: Union of k-mer hash values in all subgraphs.
     """
     # a dict mapping node to penalty for faster lookup
     node_penalty: dict[int, float] = dict(
@@ -159,9 +159,9 @@ def get_subgraphs(
         logger.info(f' - Found {len(subgraphs)} low-penalty subgraphs')
     else:
         log_and_raise(
-            RuntimeError, 
+            RuntimeError,
             ('No low-penalty subgraph was found. '
-            'Try decrease --stringency, or increase --penalty-th (penalty threshold, check log for the calculated value).')
+            'Try decrease --stringency, or increase --penalty-th (penalty threshold, check log for the calculated value)')
         )
 
     # due to the greedy nature of node expansion, subgraphs created first are usually larger
diff --git a/src/seqwin/kmers.py b/src/seqwin/kmers.py
index c7193bc..bf1f4aa 100644
--- a/src/seqwin/kmers.py
+++ b/src/seqwin/kmers.py
@@ -2,7 +2,7 @@
 K-mer Graph
 ===========
 
-A core module of Seqwin. Create a k-mer graph from k-mers of all input genome assemblies. 
+A core module of Seqwin. Create a k-mer graph from k-mers of all input genome assemblies.
 
 Dependencies:
 -------------
@@ -45,21 +45,20 @@
 
 class KmerGraph(object):
     """
-    1. Create a weighted, undirected k-mer graph, and calculate node penalty scores. 
-    2. Extract low-penalty subgraphs from the k-mer graph with `self.filter()`. 
+    1. Create a weighted, undirected k-mer graph, and calculate node penalty scores.
+    2. Extract low-penalty subgraphs from the k-mer graph with `self.filter()`.
 
     Attributes:
-        kmers (NDArray[np.void]): A 1-D NumPy structured array of k-mers from all assemblies, with dtype `KMER_DTYPE` defined in `seqwin.graph`. 
-            Grouped and sorted by k-mer hashes. 
-        idx (NDArray[np.uint64] | None): The original indices assigned when k-mers are generated (ordered by genomic positions). 
-            Parallel to `kmers`. 
-        nodes (NDArray[np.void]): A 1-D NumPy structured array of k-mer nodes, with dtype `NODE_DTYPE` defined in `seqwin.graph`. 
-            For each node, `kmers[node['start']:node['stop']]` is the k-mer group with `node['hash']`. 
-        edges (NDArray[np.void]): A 1-D NumPy structured array of weighted, undirected edges, with dtype `EDGE_DTYPE` defined in `seqwin.graph`. 
-            Edge weight is the number of assemblies where the two k-mers are adjacent. 
-        graph (nx.Graph): The graph instance built from filtered nodes and edges. 
-        subgraphs (tuple[frozenset[np.uint64], ...] | None): Low-penalty subgraphs. Each subgraph is a set of k-mer hash values. 
-            Generated with `self.filter()`. 
+        kmers (NDArray[np.void]): A 1-D NumPy structured array of k-mers from all assemblies, grouped and sorted by k-mer hashes.
+        idx (NDArray[np.uint64] | None): The original indices assigned when k-mers are generated (ordered by genomic positions).
+            Parallel to `kmers`.
+        nodes (NDArray[np.void]): A 1-D NumPy structured array of k-mer nodes.
+            For each node, `kmers[node['start']:node['stop']]` is the k-mer group with `node['hash']`.
+        edges (NDArray[np.void]): A 1-D NumPy structured array of weighted, undirected edges.
+            Edge weight is the number of assemblies where the two k-mers are adjacent.
+        graph (nx.Graph): The graph instance built from filtered nodes and edges.
+        subgraphs (tuple[frozenset[np.uint64], ...] | None): Low-penalty subgraphs. Each subgraph is a set of k-mer hash values.
+            Generated with `self.filter()`.
     """
     __slots__ = (
         'kmers', 'idx', 'nodes', 'edges', 'graph', 'subgraphs', '_filtered_flag'
@@ -74,32 +73,32 @@ class KmerGraph(object):
 
     def __init__(self, assemblies: Assemblies, kmerlen: int, windowsize: int, n_cpu: int) -> None:
         """
-        1. Generate minimizer sketches and weighted edges. 
-        2. Generate k-mer nodes and calculate node penalty scores. 
+        1. Generate minimizer sketches and weighted edges.
+        2. Generate k-mer nodes and calculate node penalty scores.
 
         Args:
-            assemblies (Assemblies): See `Assemblies` in `assemblies.py`. 
-            kmerlen (int): See `Config` in `config.py`. 
-            windowsize (int): See `Config` in `config.py`. 
-            n_cpu (int): See `Config` in `config.py`. 
+            assemblies (Assemblies): See `Assemblies` in `assemblies.py`.
+            kmerlen (int): See `Config` in `config.py`.
+            windowsize (int): See `Config` in `config.py`.
+            n_cpu (int): See `Config` in `config.py`.
         """
         n_assemblies = len(assemblies)
         logger.info(f'Building minimizer graph from {n_assemblies} assemblies...')
         tik = time()
 
         kmers, idx, nodes, edges, record_ids = build(
-            assemblies.path, 
-            kmerlen, 
-            windowsize, 
-            assemblies.index, 
-            assemblies.is_target, 
-            n_cpu=n_cpu, 
+            assemblies.path,
+            kmerlen,
+            windowsize,
+            assemblies.index,
+            assemblies.is_target,
+            n_cpu=n_cpu,
         )
         # calculate penalty for each node
         n_tar = sum(assemblies.is_target)
         n_neg = n_assemblies - n_tar
         nodes['penalty'] = _frac_to_penalty(
-            nodes['n_tar'] / n_tar, 
+            nodes['n_tar'] / n_tar,
             nodes['n_neg'] / n_neg
         )
         assemblies.record_ids = record_ids
@@ -121,17 +120,17 @@ def __init__(self, assemblies: Assemblies, kmerlen: int, windowsize: int, n_cpu:
     def filter(
         self, penalty_th: float, edge_weight_th: float, min_nodes: int, max_nodes: int | None, rng: Random
     ) -> None:
-        """Filter graph edges and nodes. 
-        1. Remove low-weight edges and isolated nodes. 
-        2. Create the graph instance and extract low-penalty subgraphs. 
-        3. Remove k-mers not included in any of the subgraphs. 
+        """Filter graph edges and nodes.
+        1. Remove low-weight edges and isolated nodes.
+        2. Create the graph instance and extract low-penalty subgraphs.
+        3. Remove k-mers not included in any of the subgraphs.
 
         Args:
-            penalty_th (float): See `RunState` in `config.py`. 
-            edge_weight_th (float): See `RunState` in `config.py`. 
-            min_nodes (int): See `RunState` in `config.py`. 
-            max_nodes (int | None): See `RunState` in `config.py`. 
-            rng (random.Random): See `RunState` in `config.py`. 
+            penalty_th (float): See `RunState` in `config.py`.
+            edge_weight_th (float): See `RunState` in `config.py`.
+            min_nodes (int): See `RunState` in `config.py`.
+            max_nodes (int | None): See `RunState` in `config.py`.
+            rng (random.Random): See `RunState` in `config.py`.
         """
         kmers = self.kmers
         idx = self.idx
@@ -173,18 +172,18 @@ def filter(
 
     @staticmethod
     def __filter_graph(nodes: NDArray, edges: NDArray, edge_weight_th: float) -> tuple[NDArray, NDArray, nx.Graph]:
-        """Remove low-weight edges and isolated nodes, and create the graph instance. 
-        
+        """Remove low-weight edges and isolated nodes, and create the graph instance.
+
         Args:
-            nodes (NDArray): See `KmerGraph.nodes`. 
-            edges (NDArray): See `KmerGraph.edges`. 
-            edge_weight_th (float): See `RunState` in `config.py`. 
+            nodes (NDArray): See `KmerGraph.nodes`.
+            edges (NDArray): See `KmerGraph.edges`.
+            edge_weight_th (float): See `RunState` in `config.py`.
 
         Returns:
             tuple: A tuple containing
-                1. NDArray: Filtered nodes. 
-                2. NDArray: Filtered edges. 
-                3. nx.Graph: See `KmerGraph.graph`. 
+                1. NDArray: Filtered nodes.
+                2. NDArray: Filtered edges.
+                3. nx.Graph: See `KmerGraph.graph`.
         """
         logger.info(' - Filtering graph edges and nodes...')
         n_nodes, n_edges = len(nodes), len(edges)
@@ -206,8 +205,8 @@ def __filter_graph(nodes: NDArray, edges: NDArray, edge_weight_th: float) -> tup
         graph = nx.Graph()
         graph.add_weighted_edges_from(edge_values, weight=EDGE_W)
         nx.set_node_attributes(
-            graph, 
-            values=dict(zip(nodes['hash'], nodes['penalty'])), 
+            graph,
+            values=dict(zip(nodes['hash'], nodes['penalty'])),
             name=NODE_P
         )
 
@@ -215,20 +214,20 @@ def __filter_graph(nodes: NDArray, edges: NDArray, edge_weight_th: float) -> tup
 
 
 def _expected_frac(jaccard_mtx: NDArray) -> np.floating:
-    """Calculate the expected fraction from a matrix of pairwise Jaccard indices. 
-    Here fraction means: for a k-mer `h` in a group of genomes (k-mer sets), the fraction of sets in 
-    a second group that also contain `h`. `jaccard_mtx` is the pairwise Jaccard indices between sets 
-    in the two groups. Note that this could be a self comparison (two groups are the same). 
-    - `E(frac) = mean(2J / (1+J))`, where `J` is the Jaccard matrix. 
-    - This expectation ≥ mean of the Jaccard matrix, (`2J / (1+J) ≥ J, 0 ≤ J ≤ 1`). 
+    """Calculate the expected fraction from a matrix of pairwise Jaccard indices.
+    Here fraction means: for a k-mer `h` in a group of genomes (k-mer sets), the fraction of sets in
+    a second group that also contain `h`. `jaccard_mtx` is the pairwise Jaccard indices between sets
+    in the two groups. Note that this could be a self comparison (two groups are the same).
+    - `E(frac) = mean(2J / (1+J))`, where `J` is the Jaccard matrix.
+    - This expectation ≥ mean of the Jaccard matrix, (`2J / (1+J) ≥ J, 0 ≤ J ≤ 1`).
     """
     return np.mean(2 * jaccard_mtx / (1 + jaccard_mtx))
 
 
 def _frac_to_penalty(frac_tar: NDArray | float, frac_neg: NDArray | float) -> NDArray | float:
-    """The penalty formula (Euclidean / L2 norm of the two fractions). 
-    - When `frac_tar` and `frac_neg` are expectations, this formula doesn't give the expected penalty, since it's non-linear. 
-    - However, since it is convex, the returned value ≤ the true expectation (Jensen’s inequality). 
+    """The penalty formula (Euclidean / L2 norm of the two fractions).
+    - When `frac_tar` and `frac_neg` are expectations, this formula doesn't give the expected penalty, since it's non-linear.
+    - However, since it is convex, the returned value ≤ the true expectation (Jensen’s inequality).
     """
     return ((1 - frac_tar)**2 + frac_neg**2)**0.5
 
@@ -236,17 +235,17 @@ def _frac_to_penalty(frac_tar: NDArray | float, frac_neg: NDArray | float) -> ND
 def get_kmers(
     assemblies: Assemblies, config: Config, state: RunState
 ) -> tuple[KmerGraph, NDArray | None]:
-    """Create a `KmerGraph` instance, calculate filtering thresholds and run the `filter()` method to generate low-penalty subgraphs. 
+    """Create a `KmerGraph` instance, calculate filtering thresholds and run the `filter()` method to generate low-penalty subgraphs.
 
     Args:
-        assemblies (Assemblies): See `Assemblies` in `assemblies.py`. 
-        config (Config): See `Config` in `config.py`. 
-        state (RunState): See `RunState` in `config.py`. 
+        assemblies (Assemblies): See `Assemblies` in `assemblies.py`.
+        config (Config): See `Config` in `config.py`.
+        state (RunState): See `RunState` in `config.py`.
 
     Returns:
         tuple: A tuple containing
-            1. KmerGraph: The KmerGraph instance. 
-            2. NDArray | None: A matrix of Jaccard indices of all assembly pairs. 
+            1. KmerGraph: The KmerGraph instance.
+            2. NDArray | None: A matrix of Jaccard indices of all assembly pairs.
     """
     overwrite = config.overwrite
     kmerlen = config.kmerlen
@@ -275,11 +274,11 @@ def get_kmers(
         # skip kmers.filter(), debug only
         kmers_path = working_dir / 'kmers.npz'
         np.savez(
-            kmers_path, 
-            kmers=kmers.kmers, 
-            idx=kmers.idx, 
-            nodes=kmers.nodes, 
-            edges=kmers.edges, 
+            kmers_path,
+            kmers=kmers.kmers,
+            idx=kmers.idx,
+            nodes=kmers.nodes,
+            edges=kmers.edges,
             allow_pickle=False
         )
         log_and_raise(SystemExit, f'Filtering is turned off. Raw k-mers are saved to {kmers_path}')
@@ -293,10 +292,10 @@ def get_kmers(
         # we only care about the presence & absence of k-mers in target assemblies
         if run_mash and HAS_MASH:
             jaccard = assemblies.mash(
-                kmerlen=kmerlen, 
-                sketchsize=sketchsize, 
-                out_path=working_dir / WORKINGDIR.mash, 
-                overwrite=overwrite, 
+                kmerlen=kmerlen,
+                sketchsize=sketchsize,
+                out_path=working_dir / WORKINGDIR.mash,
+                overwrite=overwrite,
                 n_cpu=n_cpu
             )
             e_absence_tar = 1 - _expected_frac(jaccard[:n_tar, :n_tar])
diff --git a/src/seqwin/markers.py b/src/seqwin/markers.py
index 448b8b5..bf8b837 100644
--- a/src/seqwin/markers.py
+++ b/src/seqwin/markers.py
@@ -2,8 +2,8 @@
 Markers
 =======
 
-A core module of Seqwin. Extract candidate markers (signatures) from subgraphs of a k-mer graph 
-(`kmers.KmerGraph.subgraphs`). 
+A core module of Seqwin. Extract candidate markers (signatures) from subgraphs of a k-mer graph
+(`kmers.KmerGraph.subgraphs`).
 
 Dependencies:
 -------------
@@ -49,7 +49,7 @@
 from .kmers import KmerGraph
 from .ncbi import blast
 from .graph import OrderedKmers
-from .utils import print_time_delta, log_and_raise, file_to_write, mp_wrapper, most_common, most_common_weighted
+from .utils import print_time_delta, log_and_raise, file_to_write, mp_wrapper
 from .config import Config, RunState, HAS_BLAST, WORKINGDIR, BLASTCONFIG, CONSEC_KMER_TH, LEN_TH_MUL
 
 # Set ConnectedKmers.is_bad as True if any of these warnings is present
@@ -65,18 +65,18 @@
 @dataclass(slots=True, frozen=True)
 class MarkerMetrics:
     """
-    Metrics of a marker, calculated from its BLAST alignments against target / non-target assemblies. 
-    Metrics default to None if BLAST is not run. 
+    Metrics of a marker, calculated from its BLAST alignments against target / non-target assemblies.
+    Metrics default to None if BLAST is not run.
 
     Attributes:
-        conservation (float | None): Average fraction of identical bases between the marker and target assemblies. 
-        f_tar_hits (float | None): Fraction of target assemblies with a BLAST hit. 
-        divergence (float | None): Average fraction of mismatches and gaps between the marker and non-target assemblies. 
-        f_neg_hits (float | None): Fraction of non-target assemblies with a BLAST hit. 
-        avg_repeats_tar (float | None): Average number of repeats of this marker in target assemblies. 
-        avg_pident_tar (float | None): Average percentage of identical bases of all repeats in target assemblies. 
-        avg_repeats_neg (float | None): Average number of repeats of this marker in non-target assemblies. 
-        avg_pident_neg (float | None): Average percentage of identical bases of all repeats in non-target assemblies. 
+        conservation (float | None): Average fraction of identical bases between the marker and target assemblies.
+        f_tar_hits (float | None): Fraction of target assemblies with a BLAST hit.
+        divergence (float | None): Average fraction of mismatches and gaps between the marker and non-target assemblies.
+        f_neg_hits (float | None): Fraction of non-target assemblies with a BLAST hit.
+        avg_repeats_tar (float | None): Average number of repeats of this marker in target assemblies.
+        avg_pident_tar (float | None): Average percentage of identical bases of all repeats in target assemblies.
+        avg_repeats_neg (float | None): Average number of repeats of this marker in non-target assemblies.
+        avg_pident_neg (float | None): Average percentage of identical bases of all repeats in non-target assemblies.
     """
     conservation: float | None = None
     f_tar_hits: float | None = None
@@ -94,28 +94,28 @@ class MarkerMetrics:
 
 
 class ConnectedKmers(object):
-    """The candidate marker Class, created from a low-penalty subgraph of the k-mer graph `KmerGraph.graph`. 
+    """The candidate marker Class, created from a low-penalty subgraph of the k-mer graph `KmerGraph.graph`.
 
     Attributes:
-        graph (nx.Graph): A low-penalty subgraph of the k-mer graph `KmerGraph.graph`. 
-        kmers (pd.DataFrame): K-mers of each node in the subgraph, from all assemblies. 
-            It's a subset of `KmerGraph.kmers`, with index inherited. 
-            K-mers with adjacent indices are also adjacent in the assembly sequence. 
-        loc (pd.DataFrame): Location of the subgraph in each assembly. 
-            Columns: ['assembly_idx', 'record_idx', 'start', 'stop', 'n_kmers', 
-            'kmers', 'is_target', 'n_repeats', 'len', 'seq']. 
-        path (OrderedKmers | None): K-mer ordering in the graph. None if the graph is not linear. 
-        rep (pd.Series): Representative sequence of the subgraph (a certain row in `loc`). 
-        len (int): Length of the representative sequence. 
-        n_rep (int): Number of assemblies having the same k-mer order as the representative. 
-        blast (pd.DataFrame | None): The best BLAST hit of the representative sequence in each assembly. 
-        metrics (MarkerMetrics): Metrics calculated with BLAST. 
-        rep_ratio (float | None): Fraction of target assemblies that have the same k-mer ordering as the representative. 
-        warnings (set): Undesirable features of the k-mer ordering. 
-        is_bad (bool): Set as True if `warnings` has anything listed in `_BAD_WARNINGS`. 
+        graph (nx.Graph): A low-penalty subgraph of the k-mer graph `KmerGraph.graph`.
+        kmers (pd.DataFrame): K-mers of each node in the subgraph, from all assemblies.
+            It's a subset of `KmerGraph.kmers`, with index inherited.
+            K-mers with adjacent indices are also adjacent in the assembly sequence.
+        loc (pd.DataFrame): Location of the subgraph in each assembly.
+            Columns: ['assembly_idx', 'record_idx', 'start', 'stop', 'n_kmers',
+            'kmers', 'is_target', 'n_repeats', 'len', 'seq'].
+        path (OrderedKmers | None): K-mer ordering in the graph. None if the graph is not linear.
+        rep (pd.Series): Representative sequence of the subgraph (a certain row in `loc`).
+        len (int): Length of the representative sequence.
+        n_rep (int): Number of assemblies having the same k-mer order as the representative.
+        blast (pd.DataFrame | None): The best BLAST hit of the representative sequence in each assembly.
+        metrics (MarkerMetrics): Metrics calculated with BLAST.
+        rep_ratio (float | None): Fraction of target assemblies that have the same k-mer ordering as the representative.
+        warnings (set): Undesirable features of the k-mer ordering.
+        is_bad (bool): Set as True if `warnings` has anything listed in `_BAD_WARNINGS`.
     """
     __slots__ = (
-        'graph', 'kmers', 'loc', 'path', 'rep', 'len', 'n_rep', 'blast', 
+        'graph', 'kmers', 'loc', 'path', 'rep', 'len', 'n_rep', 'blast',
         'metrics', 'rep_ratio', 'warnings', 'is_bad'
     )
     graph: nx.Graph
@@ -132,17 +132,17 @@ class ConnectedKmers(object):
     is_bad: bool
 
     def __init__(self, graph: nx.Graph, kmers: pd.DataFrame, kmerlen: int) -> None:
-        """Given a subgraph of the k-mer graph, 
-        1. Determine the boundary of the subgraph in each assembly. 
-        2. Determine the representative k-mer order. 
-        3. (deprecated) Determine the orientation (strand, +/-) of the subgraph in each assembly. 
+        """Given a subgraph of the k-mer graph,
+        1. Determine the boundary of the subgraph in each assembly.
+        2. Determine the representative k-mer order.
+        3. (deprecated) Determine the orientation (strand, +/-) of the subgraph in each assembly.
 
         Args:
-            graph (nx.Graph): A connected low-penalty subgraph of the k-mer graph `KmerGraph.graph`. 
-            kmers (pd.DataFrame): K-mers of each node in the subgraph, from all assemblies. 
-                It's a subset of `KmerGraph.kmers`, with index inherited. 
-                K-mers with adjacent indices are also adjacent in the assembly sequence. 
-            kmerlen (int): See `Config` in `config.py`. 
+            graph (nx.Graph): A connected low-penalty subgraph of the k-mer graph `KmerGraph.graph`.
+            kmers (pd.DataFrame): K-mers of each node in the subgraph, from all assemblies.
+                It's a subset of `KmerGraph.kmers`, with index inherited.
+                K-mers with adjacent indices are also adjacent in the assembly sequence.
+            kmerlen (int): See `Config` in `config.py`.
         """
         warnings = set() # passed to methods to add warnings in-place
 
@@ -191,17 +191,17 @@ def __init__(self, graph: nx.Graph, kmers: pd.DataFrame, kmerlen: int) -> None:
 
     @staticmethod
     def __get_loc(kmers: pd.DataFrame, kmerlen: int) -> pd.DataFrame:
-        """Determine the location / boundary of the subgraph in each assembly. 
-        1. Find consecutive k-mers for each assembly. 
-        2. Determine the boundaries (start & stop) of each group of consecutive k-mers. 
-        3. Select the largest consecutive group for each assembly. 
+        """Determine the location / boundary of the subgraph in each assembly.
+        1. Find consecutive k-mers for each assembly.
+        2. Determine the boundaries (start & stop) of each group of consecutive k-mers.
+        3. Select the largest consecutive group for each assembly.
 
         Args:
-            kmers (pd.DataFrame): See `ConnectedKmers.__init__()`. 
-            kmerlen (int): See `Config` in `config.py`. 
+            kmers (pd.DataFrame): See `ConnectedKmers.__init__()`.
+            kmerlen (int): See `Config` in `config.py`.
 
         Returns:
-            pd.DataFrame: See `ConnectedKmers.loc`. 
+            pd.DataFrame: See `ConnectedKmers.loc`.
         """
         # sort by k-mer index
         # essentially the same as sorting by k-mer position ['assembly_idx', 'record_idx', 'pos']
@@ -226,10 +226,10 @@ def __get_loc(kmers: pd.DataFrame, kmerlen: int) -> pd.DataFrame:
             by=['assembly_idx', 'record_idx', 'consec_gp'], as_index=False, sort=False, observed=True
         ).agg(
             # if we only need 'pos', "groupby()['pos'].agg(['first', 'last', 'size'])" is faster
-            start=pd.NamedAgg(column='pos', aggfunc='first'), 
-            stop=pd.NamedAgg(column='pos', aggfunc='last'), 
-            n_kmers=pd.NamedAgg(column='pos', aggfunc='size'), 
-            kmers=pd.NamedAgg(column='hash', aggfunc=tuple), 
+            start=pd.NamedAgg(column='pos', aggfunc='first'),
+            stop=pd.NamedAgg(column='pos', aggfunc='last'),
+            n_kmers=pd.NamedAgg(column='pos', aggfunc='size'),
+            kmers=pd.NamedAgg(column='hash', aggfunc=tuple),
             # kmer_strand=pd.NamedAgg(column='strand', aggfunc='sum'), # much faster than "lambda x: ''.join(x)"
             is_target=pd.NamedAgg(column='is_target', aggfunc='first'), # should be the same within each group
         )
@@ -255,18 +255,18 @@ def __get_loc(kmers: pd.DataFrame, kmerlen: int) -> pd.DataFrame:
 
     @staticmethod
     def __get_rep_order(loc: pd.DataFrame, warnings: set) -> tuple[OrderedKmers, int]:
-        """Determine the representative k-mer order and the number of target assemblies having it. 
-        1. Find the most common canonical k-mer ordering in target assemblies, weighted by the number of k-mers. 
-        2. Sanity check. 
+        """Determine the representative k-mer order and the number of target assemblies having it.
+        1. Find the most common canonical k-mer ordering in target assemblies, weighted by the number of k-mers.
+        2. Sanity check.
 
         Args:
-            loc (pd.DataFrame): See `ConnectedKmers.loc`. 
-            warnings (set): See `ConnectedKmers.warnings`. 
+            loc (pd.DataFrame): See `ConnectedKmers.loc`.
+            warnings (set): See `ConnectedKmers.warnings`.
 
         Returns:
             tuple: A tuple containing
-                1. OrderedKmers: The representative k-mer order. 
-                2. int: See `ConnectedKmers.n_rep`. 
+                1. OrderedKmers: The representative k-mer order.
+                2. int: See `ConnectedKmers.n_rep`.
         """
         # count the number of each unique k-mer ordering in target assemblies
         tar_kmers = loc[loc['is_target'] == True]['kmers']
@@ -281,12 +281,12 @@ def __get_rep_order(loc: pd.DataFrame, warnings: set) -> tuple[OrderedKmers, int
 
         # get the most common canonical ordering, weighted by the number of k-mers
         rep_canonical = max(
-            c_canonical, 
+            c_canonical,
             key=lambda k: len(k)*c_canonical[k]
         )
         # get the most common orientation
         rep_order = OrderedKmers(max(
-            (rep_canonical, rep_canonical[::-1]), 
+            (rep_canonical, rep_canonical[::-1]),
             key=lambda k: c[k] # if k does not exist in tar_kmers, c will return 0 (e.g., only one orientation exists)
         ))
 
@@ -300,17 +300,17 @@ def __get_rep_order(loc: pd.DataFrame, warnings: set) -> tuple[OrderedKmers, int
 
     @staticmethod
     def __get_graph_order(graph: nx.Graph, rep_order: OrderedKmers, warnings: set) -> OrderedKmers | None:
-        """Determine k-mer order in the subgraph. 
-        1. Check if the subgraph is linear. Return None if not linear. 
-        2. If linear, determine its k-mer ordering and check if it has the same orientation as `rep_order`. 
+        """Determine k-mer order in the subgraph.
+        1. Check if the subgraph is linear. Return None if not linear.
+        2. If linear, determine its k-mer ordering and check if it has the same orientation as `rep_order`.
 
         Args:
-            graph (nx.Graph): See `ConnectedKmers.__init__()`. 
-            rep_order (OrderedKmers): See `ConnectedKmers.__get_rep_order()`. 
-            warnings (set): See `ConnectedKmers.warnings`. 
+            graph (nx.Graph): See `ConnectedKmers.__init__()`.
+            rep_order (OrderedKmers): See `ConnectedKmers.__get_rep_order()`.
+            warnings (set): See `ConnectedKmers.warnings`.
 
         Returns:
-            OrderedKmers | None: If the graph is linear, return the k-mer order in the graph; else return None. 
+            OrderedKmers | None: If the graph is linear, return the k-mer order in the graph; else return None.
         """
         # check linearity
         leaf_nodes = tuple(node for node in graph if graph.degree[node] == 1)
@@ -352,144 +352,11 @@ def __get_graph_order(graph: nx.Graph, rep_order: OrderedKmers, warnings: set) -
 
         return graph_order
 
-    @staticmethod
-    def __get_strand(loc: pd.DataFrame, ref_order: OrderedKmers, warnings: set) -> pd.DataFrame:
-        """Determine the orientation of the subgraph in each assembly, based on k-mer ordering 
-        ('+': forward, '-': reverse, '?': undetermined, 'u': single k-mer). 
-
-        Args:
-            loc (pd.DataFrame): See `ConnectedKmers.loc`. 
-            ref_order (OrderedKmers): The reference k-mer ordering (representative or graph). 
-            warnings (set): See `ConnectedKmers.warnings`. 
-
-        Returns:
-            pd.DataFrame: See `ConnectedKmers.loc` (with updated 'strand' column). 
-        """
-        # initialize
-        loc['strand'] = '?'
-
-        if ref_order != ref_order.rev:
-            # forward and reverse ordering is not the same
-            loc['strand'] = loc['kmers'].apply(ref_order.which_strand)
-            warnings.update(ref_order.warning)
-
-            # handle sequences with only one shared k-mer with ref_order (strand as 'u'), deprecated
-            # not needed if using graph_order as ref, since these are sequences with only one k-mer
-            # loc = ConnectedKmers.__get_strand_single(loc, ref_order, warnings)
-        else:
-            # k-mer ordering is reversible
-            warnings.add('rev')
-            # determine strand based on 'kmer_strand', deprecated
-            # loc = ConnectedKmers.__get_strand_ks(loc)
-
-        return loc
-
-    @staticmethod
-    def __get_strand_single(loc: pd.DataFrame, ref_order: OrderedKmers, warnings: set) -> pd.DataFrame:
-        """Handle rows in `ConnectedKmers.loc` that have only one shared k-mer with `ref_order`. 
-
-        Args:
-            loc (pd.DataFrame): See `ConnectedKmers.loc`. 
-            ref_order (OrderedKmers): The reference k-mer ordering (representative or graph). 
-            warnings (set): See `ConnectedKmers.warnings`. 
-
-        Returns:
-            pd.DataFrame: See `ConnectedKmers.loc` (with updated 'strand' column). 
-        """
-        is_single = (loc['strand'] == 'u')
-        if not is_single.any():
-            return loc
-
-        # get kmer_strand of ref_order
-        # it's likely that they all have the same kmer_strand, but just in case we find the most common one
-        loc_ref = loc[loc['kmers'] == ref_order]
-        if len(loc_ref) > 0:
-            ref_strand: str = most_common(loc_ref['kmer_strand'])
-        else:
-            # ref_order is not found in loc
-            # might happen when the most common order and graph_order is not the same
-            return loc
-        strand_map = {kmer: strand for kmer, strand in zip(ref_order, ref_strand)}
-
-        # determine strand by k-mer strand
-        def which_strand(row: pd.Series) -> str:
-            """To be applied to all rows of loc[is_single]. 
-            """
-            # here row['kmers'] can have one k-mer, or multiple k-mers but only one of them is found in ref_order
-            for kmer, strand in zip(row['kmers'], row['kmer_strand']):
-                try:
-                    if strand == strand_map[kmer]:
-                        return '+'
-                    else:
-                        return '-'
-                except KeyError:
-                    # the current k-mer is not included in the most common ordering
-                    continue
-            # no shared k-mer with mc_order, which should not happen since this is handled in OrderedKmers.which_strand()
-            warnings.add('single_?')
-            return '?'
-        loc.loc[is_single, 'strand'] = loc[is_single].apply(which_strand, axis=1)
-        return loc
-
-    def __get_strand_ks(loc: pd.DataFrame, warnings: set) -> pd.DataFrame:
-        """Determine sequence strand based on 'kmer_strand'. Only used when strand cannot be determined by k-mer ordering
-        (e.g., forward and reverse ordering are the same (ABA), or mc_order has only one k-mer). 
-        NOTE: forward / reverse strand might have the same k-mer strand ordering (e.g., '++--'). 
-
-        Args:
-            loc (pd.DataFrame): See `ConnectedKmers.loc`. 
-            warnings (set): See `ConnectedKmers.warnings`. 
-
-        Returns:
-            pd.DataFrame: See `ConnectedKmers.loc` (with updated 'strand' column). 
-        """
-        # get the most common 'kmer_strand', weighted by the number of k-mers
-        mc_strand: str = most_common_weighted(loc['kmer_strand'])
-        # reverse complement of mc_strand
-        mc_strand_rc = mc_strand.translate(_KMER_STRAND_COMP)[::-1]
-
-        def which_strand(kmer_strand: str) -> str:
-            if kmer_strand == mc_strand:
-                return '+'
-            elif kmer_strand == mc_strand_rc:
-                return '-'
-            elif (kmer_strand in mc_strand) and (kmer_strand not in mc_strand_rc):
-                return '+'
-            elif (kmer_strand in mc_strand_rc) and (kmer_strand not in mc_strand):
-                return '-'
-            else: # undetermined strand
-                warnings.add('rev_?')
-                return '?'
-        loc['strand'] = loc['kmer_strand'].apply(which_strand)
-        return loc
-
-    @staticmethod
-    def __filter(loc: pd.DataFrame) -> pd.DataFrame:
-        """Remove abnormal rows in loc. 
-
-        Args:
-            loc (pd.DataFrame): See `ConnectedKmers.loc`. 
-
-        Returns:
-            pd.DataFrame: See `ConnectedKmers.loc` (with rows removed). 
-        """
-        # remove sequences with
-        # 1. only one k-mer
-        loc = loc[loc['n_kmers'] > 1]
-        # 2. undetermined strand
-        loc = loc[(loc['strand'] == '+') | (loc['strand'] == '-')]
-        # 3. abnormal length (should be done last)
-        # a possible scenario resulting in a super long sequence is a long run of 'N's in the original sequence
-        # Indexlr will skip those 'N's and there will be a huge gap between the coordinates of neighboring k-mers
-        med_len = np.median(loc['len'])
-        loc = loc[loc['len'] < LEN_TH_MUL*med_len]
-        return loc
-
 
 def _create_ck(
     graph: nx.Graph, nodes: tuple[np.uint64], kmers: tuple, idx: tuple, n_tar: int, kmerlen: int
 ) -> ConnectedKmers:
-    """Create a ConnectedKmers instance by taking the outputs of `_get_create_ck_args()`. 
+    """Create a ConnectedKmers instance by taking the outputs of `_get_create_ck_args()`.
     """
     # add hash and idx to each group
     kmers_df = list()
@@ -507,15 +374,15 @@ def _create_ck(
 def _get_create_ck_args(
     kg: KmerGraph, kmerlen: int, n_tar: int
 ) -> Generator[tuple[nx.Graph, pd.DataFrame, int], None, None]:
-    """Generate input arguments for `_create_ck()`. 
+    """Generate input arguments for `_create_ck()`.
 
     Args:
-        kg (KmerGraph): See `KmerGraph` in `kmers.py`. 
-        kmerlen (int): See `Config` in `config.py`. 
-        n_tar (int): See `RunState` in `config.py`. 
+        kg (KmerGraph): See `KmerGraph` in `kmers.py`.
+        kmerlen (int): See `Config` in `config.py`.
+        n_tar (int): See `RunState` in `config.py`.
 
     Yields:
-        tuple: Input arguments of `_create_ck()`. 
+        tuple: Input arguments of `_create_ck()`.
     """
     kmers = kg.kmers
     idx = kg.idx
@@ -550,22 +417,22 @@ def _get_create_ck_args(
 def _fetch_cks_seq(
     all_cks: list[ConnectedKmers], assemblies: Assemblies, rep_only: bool, n_cpu: int
 ) -> list[str] | None:
-    """Fetch the actual sequences for a list of ConnectedKmers instances. 
+    """Fetch the actual sequences for a list of ConnectedKmers instances.
 
     Args:
-        all_cks (list[ConnectedKmers]): See `_get_cks()`. 
-        assemblies (Assemblies): See `Assemblies` in `assemblies.py`. 
-        rep_only (bool): If True, only fetch the representative of each ConnectedKmers instance (fewer assemblies to be loaded); 
-            else fetch all sequences. 
-        n_cpu (int): See `Config` in `config.py`. 
+        all_cks (list[ConnectedKmers]): See `_get_cks()`.
+        assemblies (Assemblies): See `Assemblies` in `assemblies.py`.
+        rep_only (bool): If True, only fetch the representative of each ConnectedKmers instance (fewer assemblies to be loaded);
+            else fetch all sequences.
+        n_cpu (int): See `Config` in `config.py`.
 
     Returns:
-        list[str] | None: If `rep_only=True`, return a list of representative sequences; else return None. 
+        list[str] | None: If `rep_only=True`, return a list of representative sequences; else return None.
     """
     if rep_only:
         # concat all ConnectedKmers.rep and transpose
         df_loc = pd.concat(
-            (ck.rep for ck in all_cks), 
+            (ck.rep for ck in all_cks),
             axis=1, ignore_index=True
         ).transpose()
     else:
@@ -573,7 +440,7 @@ def _fetch_cks_seq(
         ck_idx = range(len(all_cks))
         # ck.loc.index is already sorted
         df_loc = pd.concat(
-            (ck.loc for ck in all_cks), 
+            (ck.loc for ck in all_cks),
             ignore_index=False, keys=ck_idx
         )
 
@@ -597,22 +464,22 @@ def _get_cks(
     kmers: KmerGraph, assemblies: Assemblies, kmerlen: int, min_len: int, n_tar: int, n_cpu: int
 ) -> tuple[list[ConnectedKmers], list[str]]:
     """
-    1. Create a ConnectedKmers instance for each low-penalty subgraph of the k-mer graph (`KmerGraph.subgraphs`). 
-    2. Remove instances that are shorter than min_len or have defects (`ConnectedKmers.is_bad`). 
-    3. Fetch the representative sequence for each remaining instance. 
+    1. Create a ConnectedKmers instance for each low-penalty subgraph of the k-mer graph (`KmerGraph.subgraphs`).
+    2. Remove instances that are shorter than min_len or have defects (`ConnectedKmers.is_bad`).
+    3. Fetch the representative sequence for each remaining instance.
 
     Args:
-        kmers (KmerGraph): See `KmerGraph` in `kmers.py`. 
-        assemblies (Assemblies): See `Assemblies` in `assemblies.py`. 
-        kmerlen (int): See `Config` in `config.py`. 
-        min_len (int): See `Config` in `config.py`. 
-        n_tar (int): See `RunState` in `config.py`. 
-        n_cpu (int): See `Config` in `config.py`. 
+        kmers (KmerGraph): See `KmerGraph` in `kmers.py`.
+        assemblies (Assemblies): See `Assemblies` in `assemblies.py`.
+        kmerlen (int): See `Config` in `config.py`.
+        min_len (int): See `Config` in `config.py`.
+        n_tar (int): See `RunState` in `config.py`.
+        n_cpu (int): See `Config` in `config.py`.
 
     Returns:
         tuple: A tuple containing
-            1. list[ConnectedKmers]: Candidate markers as ConnectedKmers instances. 
-            2. list[str]: Sequences of each marker. 
+            1. list[ConnectedKmers]: Candidate markers as ConnectedKmers instances.
+            2. list[str]: Sequences of each marker.
     """
     logger.info('Finding a representative for each low-penalty subgraph...')
     tik = time()
@@ -620,14 +487,14 @@ def _get_cks(
     # create a ConnectedKmers instance for each subgraph
     logger.info(' - Processing each subgraph...')
     all_cks: list[ConnectedKmers] = mp_wrapper(
-        _create_ck, 
-        _get_create_ck_args(kmers, kmerlen, n_tar), 
+        _create_ck,
+        _get_create_ck_args(kmers, kmerlen, n_tar),
         n_cpu=n_cpu, n_jobs=len(kmers.subgraphs)
     )
 
     # get candidate ConnectedKmers instances
     all_cks = list(
-        ck for ck in all_cks 
+        ck for ck in all_cks
         if (ck.len >= min_len) and (not ck.is_bad)
     )
     logger.info(f' - Found {len(all_cks)} candidate signatures')
@@ -644,36 +511,36 @@ def _get_cks(
 
 
 def _get_avg_ident(blast_out: pd.DataFrame, query_len: int, n: int) -> float:
-    """Given a list of BLAST hits, calculate the average sequence identity between the query and all subjects. 
-    The denominator (`n`) is the number of subject sequences that are expected to include the query sequence. 
-    Note that `n` might not be the same as `len(blast_out)`, since some subjects may have no hit. 
+    """Given a list of BLAST hits, calculate the average sequence identity between the query and all subjects.
+    The denominator (`n`) is the number of subject sequences that are expected to include the query sequence.
+    Note that `n` might not be the same as `len(blast_out)`, since some subjects may have no hit.
 
     Args:
         blast_out (pd.DataFrame): Each row should be a BLAST hit of the query, with column
-            'nident' (number of identical matches). 
-        query_len (int): Length of the query sequence. 
-        n (int): The number of subjects that are expected to include the query sequence. 
+            'nident' (number of identical matches).
+        query_len (int): Length of the query sequence.
+        n (int): The number of subjects that are expected to include the query sequence.
 
     Returns:
-        float: Conservation. 
+        float: Conservation.
     """
     return sum(blast_out['nident']) / query_len / n
 
 
 def _get_avg_dist(blast_out: pd.DataFrame, query_len: int, n: int) -> float:
-    """Given a list of BLAST hits, calculate the average distance between the query and all subjects. 
-    The denominator (`n`) is the number of subject sequences that are expected to include the query sequence. 
-    Note that `n` might not be the same as `len(blast_out)`, since some subjects may have no hit. 
+    """Given a list of BLAST hits, calculate the average distance between the query and all subjects.
+    The denominator (`n`) is the number of subject sequences that are expected to include the query sequence.
+    Note that `n` might not be the same as `len(blast_out)`, since some subjects may have no hit.
 
     Args:
-        blast_out (pd.DataFrame): Each row should be a BLAST hit of the query, with columns, 
-            1. 'mismatch': number of mismatches. 
-            2. 'gaps': total number of gaps in BOTH query and subject (might cause inaccuracy). 
-        query_len (int): Length of the query sequence. 
-        n (int): The number of subjects that are expected to include the query sequence. 
+        blast_out (pd.DataFrame): Each row should be a BLAST hit of the query, with columns,
+            1. 'mismatch': number of mismatches.
+            2. 'gaps': total number of gaps in BOTH query and subject (might cause inaccuracy).
+        query_len (int): Length of the query sequence.
+        n (int): The number of subjects that are expected to include the query sequence.
 
     Returns:
-        float: Divergence. 
+        float: Divergence.
     """
     return sum(blast_out['mismatch'] + blast_out['gaps']) / query_len / n
 
@@ -681,19 +548,19 @@ def _get_avg_dist(blast_out: pd.DataFrame, query_len: int, n: int) -> float:
 def _get_metrics(
     blast_out: pd.DataFrame, marker_len: int, n_tar: int, n_neg: int
 ) -> MarkerMetrics:
-    """Calculate the metrics of a marker based on its BLAST hits in all assemblies. 
-    - Conservation is calculated with `_get_avg_ident()` on target assemblies. 
-    - Divergence is calculated with `_get_avg_dist()` on non-target assemblies. 
+    """Calculate the metrics of a marker based on its BLAST hits in all assemblies.
+    - Conservation is calculated with `_get_avg_ident()` on target assemblies.
+    - Divergence is calculated with `_get_avg_dist()` on non-target assemblies.
 
     Args:
-        blast_out (pd.DataFrame): Each row is the best BLAST hit of the marker in an assembly. 
+        blast_out (pd.DataFrame): Each row is the best BLAST hit of the marker in an assembly.
             Required columns: ['is_target', 'nident', 'mismatch', 'gaps', 'n_hits', 'avg_nident']
-        marker_len (int): Marker length. 
-        n_tar (int): Number of target assemblies. 
-        n_neg (int): Number of non-target assemblies. 
+        marker_len (int): Marker length.
+        n_tar (int): Number of target assemblies.
+        n_neg (int): Number of non-target assemblies.
 
     Returns:
-        MarkerMetrics: Marker metrics. 
+        MarkerMetrics: Marker metrics.
     """
     if blast_out is None: # no blast hit in any assembly
         return _BASELINE_METRICS
@@ -722,19 +589,19 @@ def _get_metrics(
 def eval_markers(
     all_seqs: list[str], blastdb: Path, n_tar: int, n_neg: int, n_cpu: int=1
 ) -> tuple[list[pd.DataFrame], list[MarkerMetrics]]:
-    """BLAST check each marker (signature) sequence against all / non-target assemblies, and calculate the metrics of each marker. 
+    """BLAST check each marker (signature) sequence against all / non-target assemblies, and calculate the metrics of each marker.
 
     Args:
-        all_seqs (list[str]): A list of marker sequences. 
-        blastdb (Path): Path to a BLAST database generated by Seqwin (e.g., `seqwin-out/blastdb/`). 
-        n_tar (int): Number of target assemblies. 
-        n_neg (int): Number of non-target assemblies. 
+        all_seqs (list[str]): A list of marker sequences.
+        blastdb (Path): Path to a BLAST database generated by Seqwin (e.g., `seqwin-out/blastdb/`).
+        n_tar (int): Number of target assemblies.
+        n_neg (int): Number of non-target assemblies.
         n_cpu (int, optional): Number of threads to use. [1]
 
     Returns:
         tuple: A tuple containing
-            1. list[pd.DataFrame]: BLAST hits of each marker. 
-            2. list[MarkerMetrics]: Metrics of each marker. 
+            1. list[pd.DataFrame]: BLAST hits of each marker.
+            2. list[MarkerMetrics]: Metrics of each marker.
     """
     if blastdb.name == BLASTCONFIG.title_neg_only:
         neg_only = True
@@ -750,7 +617,7 @@ def eval_markers(
     # blast check all markers against all / non-target assemblies
     blast_out = blast(all_seqs, db=blastdb, task=BLASTCONFIG.task, columns=BLASTCONFIG.columns, n_cpu=n_cpu, batch_size=BLASTCONFIG.batch_size)
     if len(blast_out) == 0:
-        log_and_raise(RuntimeError, 'No BLAST hit found.')
+        log_and_raise(RuntimeError, 'No BLAST hit found')
     # blast_out.to_pickle('blast_out.pkl')
 
     #---------- extract BLAST hits of each marker ----------#
@@ -766,7 +633,7 @@ def eval_markers(
 
     # keep the best alignment (highest bitscore) for each assembly
     blast_out.sort_values(
-        by=['qseqid', 'assembly_idx', 'bitscore'], 
+        by=['qseqid', 'assembly_idx', 'bitscore'],
         ascending=[True, True, False], inplace=True
     )
     blast_out = blast_out.groupby(
@@ -774,7 +641,7 @@ def eval_markers(
     )
     # also keep nident of other less optimal alignments
     nident = blast_out['nident'].agg(
-        n_hits='count', 
+        n_hits='count',
         avg_nident='mean'
     )
     blast_out = blast_out.head(1)
@@ -798,9 +665,9 @@ def eval_markers(
     # calculate conservation and divergence for each marker based on its blast output
     logger.info(' - Evaluating each signature...')
     metrics_args = zip(
-        all_blast, 
-        map(len, all_seqs), 
-        repeat(n_tar, n_seqs), 
+        all_blast,
+        map(len, all_seqs),
+        repeat(n_tar, n_seqs),
         repeat(n_neg, n_seqs)
     )
     metrics = mp_wrapper(
@@ -815,18 +682,18 @@ def _eval_cks(
     all_cks: list[ConnectedKmers], all_reps: list[str], blastdb: Path, n_tar: int, n_neg: int, n_cpu: int
 ) -> None:
     """
-    1. BLAST check the representative sequence of each ConnectedKmers instance (ck), against all / non-target assemblies. 
-    2. Calculate the conservation and divergence for each ck. 
-    3. Update the attributes of each ck. 
-    4. Sort all_cks by conservation + divergence. 
+    1. BLAST check the representative sequence of each ConnectedKmers instance (ck), against all / non-target assemblies.
+    2. Calculate the conservation and divergence for each ck.
+    3. Update the attributes of each ck.
+    4. Sort all_cks by conservation + divergence.
 
     Args:
-        all_cks (list[ConnectedKmers]): ConnectedKmers instances. 
-        all_reps (list[str]): Representative sequences, in the same order as all_cks. 
-        blastdb (Path): See `RunState` in `config.py`. 
-        n_tar (int): See `RunState` in `config.py`. 
-        n_neg (int): See `RunState` in `config.py`. 
-        n_cpu (int): See `Config` in `config.py`. 
+        all_cks (list[ConnectedKmers]): ConnectedKmers instances.
+        all_reps (list[str]): Representative sequences, in the same order as all_cks.
+        blastdb (Path): See `RunState` in `config.py`.
+        n_tar (int): See `RunState` in `config.py`.
+        n_neg (int): See `RunState` in `config.py`.
+        n_cpu (int): See `Config` in `config.py`.
     """
     # run evaluation
     results = eval_markers(all_reps, blastdb, n_tar, n_neg, n_cpu)
@@ -837,7 +704,7 @@ def _eval_cks(
 
     # sort in-place
     all_cks.sort(
-        key=lambda ck: ck.metrics.conservation+ck.metrics.divergence, 
+        key=lambda ck: ck.metrics.conservation+ck.metrics.divergence,
         reverse=True
     )
 
@@ -845,16 +712,16 @@ def _eval_cks(
 def get_markers(
     kmers: KmerGraph, assemblies: Assemblies, config: Config, state: RunState
 ) -> list[ConnectedKmers]:
-    """Extract candidate markers (signatures) from a k-mer graph, and save them to files. 
+    """Extract candidate markers (signatures) from a k-mer graph, and save them to files.
 
     Args:
-        kmers (KmerGraph): See `KmerGraph` in `kmers.py`. 
-        assemblies (Assemblies): See `Assemblies` in `assemblies.py`. 
-        config (Config): See `Config` in `config.py`. 
-        state (RunState): See `RunState` in `config.py`. 
+        kmers (KmerGraph): See `KmerGraph` in `kmers.py`.
+        assemblies (Assemblies): See `Assemblies` in `assemblies.py`.
+        config (Config): See `Config` in `config.py`.
+        state (RunState): See `RunState` in `config.py`.
 
     Returns:
-        list[ConnectedKmers]: Candidate markers. 
+        list[ConnectedKmers]: Candidate markers.
     """
     overwrite = config.overwrite
     kmerlen = config.kmerlen
@@ -875,9 +742,9 @@ def get_markers(
     if run_blast and HAS_BLAST:
         logger.info('Evaluating candidate signatures with BLAST...')
         blastdb = assemblies.makeblastdb(
-            prefix=working_dir / WORKINGDIR.blast_dir, 
-            neg_only=blast_neg_only, 
-            overwrite=overwrite, 
+            prefix=working_dir / WORKINGDIR.blast_dir,
+            neg_only=blast_neg_only,
+            overwrite=overwrite,
             n_cpu=n_cpu
         )
         _eval_cks(all_cks, all_reps, blastdb, n_tar, n_neg, n_cpu)
@@ -910,7 +777,7 @@ def get_markers(
     markers_csv = working_dir / WORKINGDIR.markers_csv
     file_to_write(markers_csv, overwrite)
     pd.DataFrame(
-        csv, 
+        csv,
         columns=('fasta_header', 'length', *_METRIC_NAMES, 'rep_ratio', 'n_nodes')
     ).to_csv(markers_csv, index=False, encoding='utf-8', lineterminator='\n')
     logger.info(f'Metrics of candidate signatures saved as {markers_csv}')
diff --git a/src/seqwin/mash.py b/src/seqwin/mash.py
index 7b1d3a1..159dd52 100644
--- a/src/seqwin/mash.py
+++ b/src/seqwin/mash.py
@@ -2,7 +2,7 @@
 Mash
 ====
 
-Assembly distance estimation with `Mash <https://doi.org/10.1186/s13059-016-0997-x>`__. 
+Assembly distance estimation with `Mash <https://doi.org/10.1186/s13059-016-0997-x>`__.
 
 Dependencies:
 ------------
@@ -37,31 +37,31 @@
 
 
 def sketch(
-    assembly_path: Path | Iterable[Path], 
-    kmerlen: int = 21, 
-    sketchsize: int = 1000, 
-    out_path: Path | None = None, 
-    overwrite: bool = False, 
+    assembly_path: Path | Iterable[Path],
+    kmerlen: int = 21,
+    sketchsize: int = 1000,
+    out_path: Path | None = None,
+    overwrite: bool = False,
     n_cpu: int = 1
 ) -> Path:
-    """Create a Mash / MinHash sketch for a single assembly, or a merged sketch for multiple assemblies. 
+    """Create a Mash / MinHash sketch for a single assembly, or a merged sketch for multiple assemblies.
 
     Args:
-        assembly_path (Path | list[Path]): Path to the assembly file in FASTA format, or a list of assembly paths. 
+        assembly_path (Path | list[Path]): Path to the assembly file in FASTA format, or a list of assembly paths.
         kmerlen (int, optional): k-mer length. [21]
         sketchsize (int, optional): Sketch size (number of non-redundant k-mers). [1000]
-        out_path (Path | None, optional): Output path for the sketch file. 
+        out_path (Path | None, optional): Output path for the sketch file.
             If None, output to `assembly_path + '.msh'` or `assembly_path[0] + '.msh'`. [None]
         overwrite (bool, optional): If True, overwrite the existing Mash sketch file. [False]
         n_cpu (int, optional): Number of processes to run in parallel. [1]
 
     Returns:
-        Path: Path to the Mash sketch file (.msh). 
+        Path: Path to the Mash sketch file (.msh).
     """
     args = [
-        'mash', 'sketch', 
-        '-k', str(kmerlen), 
-        '-s', str(sketchsize), 
+        'mash', 'sketch',
+        '-k', str(kmerlen),
+        '-s', str(sketchsize),
         '-p', str(n_cpu)
     ]
 
@@ -79,7 +79,7 @@ def sketch(
         log_text = f' - Generating MinHash sketches with Mash for {len(assembly_path)} assemblies...'
         assembly_path = assembly_path[0]
     else:
-        log_and_raise(ValueError, 'Invalid assembly_path for mash_sketch, must be str or list.')
+        log_and_raise(ValueError, 'Invalid assembly_path for mash_sketch, must be str or list')
 
     # determine output path (Mash appends .msh to the output path)
     if out_path is None:
@@ -104,30 +104,30 @@ def sketch(
 
 
 def dist(
-    ref_path: Path, 
-    query_path: Path | None = None, 
+    ref_path: Path,
+    query_path: Path | None = None,
     n_cpu: int = 1
 ) -> pd.DataFrame:
-    """Run `mash dist {ref_path} {query_path}` and parse the TSV output as a pandas DataFrame. 
-    Note: high memory usage when the number of assemblies in the sketch file is large. 
-    
+    """Run `mash dist {ref_path} {query_path}` and parse the TSV output as a pandas DataFrame.
+    Note: high memory usage when the number of assemblies in the sketch file is large.
+
     Args:
-        ref_path (Path): Path to the reference sketch file. It could include multiple sketches, merged with the 'mash paste' command. 
+        ref_path (Path): Path to the reference sketch file. It could include multiple sketches, merged with the 'mash paste' command.
         query_path (Path | None, optional): Path to the query sketch. If None, query ref_path against itself. [None]
         n_cpu (int, optional): Number of processes to run in parallel. [1]
 
     Returns:
-        pd.DataFrame: Tabular output of `mash dist`, each row is an assembly pair, 
-            with columns ['ref', 'query', 'dist', 'pval', 'jaccard', 'shared', 'total']. 
+        pd.DataFrame: Tabular output of `mash dist`, each row is an assembly pair,
+            with columns ['ref', 'query', 'dist', 'pval', 'jaccard', 'shared', 'total'].
     """
     if query_path is None:
         query_path = ref_path
 
     logger.info(' - Calculating Mash distances of assembly pairs...')
     cmd_out = run_cmd(
-        'mash', 'dist', 
-        '-p', str(n_cpu), 
-        ref_path, 
+        'mash', 'dist',
+        '-p', str(n_cpu),
+        ref_path,
         query_path
     )
     # parse stdout; pd.read_csv() does auto type conversion
@@ -138,32 +138,32 @@ def dist(
 
 
 def get_jaccard(
-    ref_path: Path, 
-    query_path: Path | None = None, 
-    n_cpu: int = 1, 
+    ref_path: Path,
+    query_path: Path | None = None,
+    n_cpu: int = 1,
     bufsize: int = 1_000_000
 ) -> Generator[float, None, None]:
-    """Estimate the pairwise Jaccard index between (each assembly sketch in the query) 
-    and (each assembly sketch in the reference), with `mash dist`. Use a stream pipe to save memory. 
-    
+    """Estimate the pairwise Jaccard index between (each assembly sketch in the query)
+    and (each assembly sketch in the reference), with `mash dist`. Use a stream pipe to save memory.
+
     Args:
-        ref_path (Path): Path to the reference sketch file. It could include multiple sketches, merged with the 'mash paste' command. 
+        ref_path (Path): Path to the reference sketch file. It could include multiple sketches, merged with the 'mash paste' command.
         query_path (Path | None, optional): Path to the query sketch. If None, query ref_path against itself. [None]
         n_cpu (int, optional): Number of processes to run in parallel. [1]
         bufsize (int, optional): `bufsize` option for `subprocess.Popen`. [1_000_000]
 
     Yields:
-        float: Jaccard index of each assembly pair. 
+        float: Jaccard index of each assembly pair.
     """
     if query_path is None:
         query_path = ref_path
 
     logger.info(' - Calculating Jaccard indices of assembly pairs...')
     proc = subprocess.Popen(
-        ('mash', 'dist', '-p', str(n_cpu), ref_path, query_path), 
-        stdout=subprocess.PIPE, 
+        ('mash', 'dist', '-p', str(n_cpu), ref_path, query_path),
+        stdout=subprocess.PIPE,
         stderr=subprocess.PIPE,
-        text=True, 
+        text=True,
         bufsize=bufsize
     )
 
diff --git a/src/seqwin/ncbi.py b/src/seqwin/ncbi.py
index 19cbe66..6a6b246 100644
--- a/src/seqwin/ncbi.py
+++ b/src/seqwin/ncbi.py
@@ -2,8 +2,8 @@
 NCBI
 ====
 
-- Search and download genome assemblies with NCBI `Datasets <https://www.ncbi.nlm.nih.gov/datasets/>`__. 
-- Run NCBI `BLAST+ <https://www.ncbi.nlm.nih.gov/books/NBK131777/>`__. 
+- Search and download genome assemblies with NCBI `Datasets <https://www.ncbi.nlm.nih.gov/datasets/>`__.
+- Run NCBI `BLAST+ <https://www.ncbi.nlm.nih.gov/books/NBK131777/>`__.
 
 Dependencies:
 -------------
@@ -99,20 +99,20 @@ def _add_api_key(args: list[str | Path], api_key: str | None) -> list[str | Path
 
 
 def search_taxon(taxon: str, api_key: str | None=None) -> tuple[str | None, str | None]:
-    """Search a taxon on NCBI Taxonomy. Internet connection is needed. 
+    """Search a taxon on NCBI Taxonomy. Internet connection is needed.
 
     Args:
-        taxon (str): Name or ID of the taxon (exact match). 
+        taxon (str): Name or ID of the taxon (exact match).
         api_key (str | None, optional): NCBI API key. [None]
 
     Returns:
         tuple: A tuple containing
-            1. str | None: NCBI Taxonomy ID of the taxon. 
-            2. str | None: Current scientific name of the taxon. 
+            1. str | None: NCBI Taxonomy ID of the taxon.
+            2. str | None: Current scientific name of the taxon.
     """
     logger.info(f'Searching NCBI Taxonomy for "{taxon}"...')
     args = [
-        'datasets', 'summary', 'taxonomy', 'taxon', str(taxon), 
+        'datasets', 'summary', 'taxonomy', 'taxon', str(taxon),
         '--as-json-lines', # output as json
         '--report', 'names', # do not output tax ids of children
     ]
@@ -133,13 +133,13 @@ def search_taxon(taxon: str, api_key: str | None=None) -> tuple[str | None, str
 
 
 def get_assembly_paths(package_dir: Path) -> list[Path]:
-    """Get the file paths of all genome assemblies in a NCBI genome package. 
+    """Get the file paths of all genome assemblies in a NCBI genome package.
 
     Args:
-        package_dir (Path): Path of the genome package directory (e.g., ncbi_dataset). 
+        package_dir (Path): Path of the genome package directory (e.g., ncbi_dataset).
 
     Returns:
-        list[Path]: File paths of all genome assemblies in the package. 
+        list[Path]: File paths of all genome assemblies in the package.
     """
     # sanity check
     if not package_dir.is_dir():
@@ -159,23 +159,23 @@ def get_assembly_paths(package_dir: Path) -> list[Path]:
 
 
 def download_taxon(
-    taxon: str, 
-    prefix: Path = Path.cwd(), 
-    format: Format = Format.fasta, 
-    level: Level = Level.contig, 
-    source: Source = Source.genbank, 
-    annotated: bool = True, 
-    exclude_mag: bool = False, 
-    gzip: bool = True, 
-    api_key: str | None = None, 
-    overwrite: bool = False, 
+    taxon: str,
+    prefix: Path = Path.cwd(),
+    format: Format = Format.fasta,
+    level: Level = Level.contig,
+    source: Source = Source.genbank,
+    annotated: bool = True,
+    exclude_mag: bool = False,
+    gzip: bool = True,
+    api_key: str | None = None,
+    overwrite: bool = False,
     n_cpu: int = 1
 ) -> list[Path] | None:
-    """Download genome assemblies under a taxon from NCBI Taxonomy. Internet connection is needed. 
-    Atypical genomes and genomes from large multi-isolate projects are excluded. 
+    """Download genome assemblies under a taxon from NCBI Taxonomy. Internet connection is needed.
+    Atypical genomes and genomes from large multi-isolate projects are excluded.
 
     Args:
-        taxon (str): Name or ID of the taxon (exact match). 
+        taxon (str): Name or ID of the taxon (exact match).
         prefix (Path, optional): A directory where the data package is downloaded. [cwd]
         format (str, optional): Format of genome sequences, should be 'fasta' or 'genbank'. ['fasta']
         level (str, optional): Limit to genomes ≥ this assembly level ('contig' < 'scaffold' < 'chromosome' < 'complete'). ['contig']
@@ -188,7 +188,7 @@ def download_taxon(
         n_cpu (int, optional): Number of processes to run in parallel. [1]
 
     Returns:
-        list[Path] | None: File paths of downloaded genome assemblies. Return None if the taxon is not found. 
+        list[Path] | None: File paths of downloaded genome assemblies. Return None if the taxon is not found.
     """
     # sanity check
     if not prefix.is_dir():
@@ -203,9 +203,9 @@ def download_taxon(
             assembly_paths = get_assembly_paths(tax_dir)
         except Exception as e:
             log_and_raise(
-                RuntimeError, 
+                RuntimeError,
                 (f'Genome package might be incomplete {tax_dir}\n'
-                'Consider deleting it and try again.'), 
+                'Consider deleting it and try again'),
                 from_e=e
             )
         logger.info(f' - Found {len(assembly_paths)} genome assemblies.')
@@ -225,11 +225,11 @@ def download_taxon(
 
     # arguments for the download command
     args = [
-        'datasets', 'download', 'genome', 
-        'taxon', tax_id, 
-        '--filename', tax_zip, 
-        '--exclude-atypical', '--exclude-multi-isolate', 
-        '--no-progressbar', '--dehydrated', 
+        'datasets', 'download', 'genome',
+        'taxon', tax_id,
+        '--filename', tax_zip,
+        '--exclude-atypical', '--exclude-multi-isolate',
+        '--no-progressbar', '--dehydrated',
     ]
     if format == Format.fasta:
         args += ['--include', 'genome']
@@ -286,7 +286,7 @@ def download_taxon(
 
     # download actual sequences
     args = [
-        'datasets', 'rehydrate', '--directory', tax_dir, 
+        'datasets', 'rehydrate', '--directory', tax_dir,
         '--max-workers', str(n_cpu), '--no-progressbar'
     ]
     if gzip:
@@ -300,10 +300,10 @@ def download_taxon(
     except Exception as e:
         shutil.rmtree(tax_dir) # remove incomplete dir
         log_and_raise(
-            RuntimeError, 
+            RuntimeError,
             (f'Failed to rehydrate data package for taxon "{taxon}".\n'
             'NCBI might have blocked the request due to high usage. Try waiting for a while before running Seqwin again.\n'
-            'Add --overwrite to the command so that downloaded taxon packages can be reused.'), 
+            'Add --overwrite to the command so that downloaded taxon packages can be reused.'),
             from_e=e
         )
 
@@ -315,8 +315,8 @@ def download_taxon(
 
 
 def _get_blast_outfmt(columns: Sequence[str]) -> str:
-    """Given the columns to be included in the BLAST TSV output, get the value to be provided to 
-    the -outfmt argument of the blastn command. See `blastn -help` for more info. 
+    """Given the columns to be included in the BLAST TSV output, get the value to be provided to
+    the -outfmt argument of the blastn command. See `blastn -help` for more info.
     """
     # we need a custom output format since some info (e.g., sequences) are not included in the tsv by default
     # do not surround the return str with double quotes (" ") even if there are spaces in it, `subprocess.run()` will do it for you
@@ -325,31 +325,31 @@ def _get_blast_outfmt(columns: Sequence[str]) -> str:
 
 
 def _blast_batch(
-    seq_idx: Sequence[int], 
-    seq_list: Sequence[str], 
-    db: Path, 
-    task: str, 
-    columns: Sequence[str], 
-    outfmt: str, 
-    taxids: str | None, 
-    neg_taxids: str | None, 
+    seq_idx: Sequence[int],
+    seq_list: Sequence[str],
+    db: Path,
+    task: str,
+    columns: Sequence[str],
+    outfmt: str,
+    taxids: str | None,
+    neg_taxids: str | None,
     n_cpu: int
 ) -> pd.DataFrame:
-    """Run BLAST on a list of sequences and return a DataFrame of the tabular output. 
+    """Run BLAST on a list of sequences and return a DataFrame of the tabular output.
 
     Args:
-        seq_idx (Sequence[int]): Indices of the input sequences. 
-        seq_list (Sequence[str]): A list of sequences for BLAST. 
-        db (Path): Path to the BLAST database. 
-        task (str): Preset BLAST tasks ('blastn', 'blastn-short', 'megablast'). 
-        columns (Sequence[str]): Columns to be included in the DataFrame output. 
-        outfmt (str): Output of `_get_blast_outfmt()`. 
-        taxids (str | None): Only search for these taxonomy IDs and their descendants. 
-        neg_taxids (str | None): Do NOT search for these taxonomy IDs and their descendants. 
-        n_cpu (int): Number of processes to run in parallel. 
+        seq_idx (Sequence[int]): Indices of the input sequences.
+        seq_list (Sequence[str]): A list of sequences for BLAST.
+        db (Path): Path to the BLAST database.
+        task (str): Preset BLAST tasks ('blastn', 'blastn-short', 'megablast').
+        columns (Sequence[str]): Columns to be included in the DataFrame output.
+        outfmt (str): Output of `_get_blast_outfmt()`.
+        taxids (str | None): Only search for these taxonomy IDs and their descendants.
+        neg_taxids (str | None): Do NOT search for these taxonomy IDs and their descendants.
+        n_cpu (int): Number of processes to run in parallel.
 
     Returns:
-        pd.DataFrame: A DataFrame of the tabular output. 
+        pd.DataFrame: A DataFrame of the tabular output.
     """
     # create input fasta for blast
     blast_in = ''.join(
@@ -358,9 +358,9 @@ def _blast_batch(
 
     # prepare blastn args
     args = [
-        'blastn', 
-        '-db', db, 
-        '-task', task, 
+        'blastn',
+        '-db', db,
+        '-task', task,
         '-outfmt', outfmt, # output in tsv format with specified columns
         '-max_hsps', _MAX_HSPS, # maximum number of HSPs per subject sequence to save for each query
         '-max_target_seqs', _MAX_TARGET_SEQS, # maximum number of aligned sequences to keep
@@ -381,24 +381,24 @@ def _blast_batch(
 
 
 def blast(
-    seq_list: Sequence[str], 
-    db: Path, 
-    task: Task = Task.blastn, 
-    columns: Sequence[str] | None = None, 
-    taxids: Sequence[int] | None = None, 
-    neg_taxids: Sequence[int] | None = None, 
-    n_cpu: int = 1, 
+    seq_list: Sequence[str],
+    db: Path,
+    task: Task = Task.blastn,
+    columns: Sequence[str] | None = None,
+    taxids: Sequence[int] | None = None,
+    neg_taxids: Sequence[int] | None = None,
+    n_cpu: int = 1,
     batch_size: int = 1000
 ) -> pd.DataFrame:
-    """Run BLAST on a list of sequences and return a DataFrame of the tabular output. 
-    - -max_hsps is set to 1000. 
-    - -max_target_seqs is set to 50000. 
+    """Run BLAST on a list of sequences and return a DataFrame of the tabular output.
+    - -max_hsps is set to 1000.
+    - -max_target_seqs is set to 50000.
     - To avoid having too many query sequences in a single BLAST run, each run only takes a batch
-    of sequences with batch size determined by batch_size. 
+    of sequences with batch size determined by batch_size.
 
     Args:
-        seq_list (Sequence[str]): A list of sequences for BLAST. 
-        db (Path): Path to the BLAST database. 
+        seq_list (Sequence[str]): A list of sequences for BLAST.
+        db (Path): Path to the BLAST database.
         task (str): Preset BLAST tasks ('blastn', 'blastn-short', 'megablast'). ['blastn']
         columns (Sequence[str] | None): Columns to be included in the DataFrame output. None to use default columns. [None]
         taxids (Sequence[int] | None): Only search for these taxonomy IDs and their descendants. [None]
@@ -407,23 +407,23 @@ def blast(
         batch_size (int, optional): Number of query sequences in a single BLAST run. [1000]
 
     Returns:
-        pd.DataFrame: A DataFrame of the tabular output, with default columns, 
-            1. 'qseqid': query or source (e.g., gene) sequence id. 
-            2. 'sseqid': subject or target (e.g., reference genome) sequence id. 
-            3. 'length': alignment length. 
-            4. 'pident': percentage of identical matches (= nident / length). 
-            5. 'nident': number of identical matches (= length - mismatch - gaps). 
-            6. 'mismatch': number of mismatches. 
-            7. 'gapopen': number of gap openings. 
-            8. 'gaps': total number of gaps in both query and subject. 
-            9. 'qstart': start of alignment in query. 
-            10. 'qend': end of alignment in query. 
-            11. 'sstart': start of alignment in subject. 
-            12. 'send': end of alignment in subject. 
-            13. 'evalue': expect value. 
-            14. 'bitscore': bit score. 
-            15. 'qseq': aligned part of query sequence. 
-            16. 'sseq': aligned part of subject sequence. 
+        pd.DataFrame: A DataFrame of the tabular output, with default columns,
+            1. 'qseqid': query or source (e.g., gene) sequence id.
+            2. 'sseqid': subject or target (e.g., reference genome) sequence id.
+            3. 'length': alignment length.
+            4. 'pident': percentage of identical matches (= nident / length).
+            5. 'nident': number of identical matches (= length - mismatch - gaps).
+            6. 'mismatch': number of mismatches.
+            7. 'gapopen': number of gap openings.
+            8. 'gaps': total number of gaps in both query and subject.
+            9. 'qstart': start of alignment in query.
+            10. 'qend': end of alignment in query.
+            11. 'sstart': start of alignment in subject.
+            12. 'send': end of alignment in subject.
+            13. 'evalue': expect value.
+            14. 'bitscore': bit score.
+            15. 'qseq': aligned part of query sequence.
+            16. 'sseq': aligned part of subject sequence.
     """
     # check input
     tot_seq = len(seq_list)
@@ -451,9 +451,9 @@ def blast(
         logger.info(f' - {batch_start}/{tot_seq}')
         batch_stop = batch_start + batch_size
         blast_out.append(_blast_batch(
-            seq_idx[batch_start: batch_stop], 
-            seq_list[batch_start: batch_stop], 
-            db, task, columns, outfmt, 
+            seq_idx[batch_start: batch_stop],
+            seq_list[batch_start: batch_stop],
+            db, task, columns, outfmt,
             taxids, neg_taxids, n_cpu
         ))
         batch_start = batch_stop
diff --git a/src/seqwin/utils.py b/src/seqwin/utils.py
index c8aa310..b5d3050 100644
--- a/src/seqwin/utils.py
+++ b/src/seqwin/utils.py
@@ -6,39 +6,6 @@
 -------------
 - numpy
 - biopython (optional)
-
-Classes:
---------
-- SharedArr
-
-Functions:
-----------
-- print_time_delta
-- log_and_raise
-- overwrite_warning
-- overwrite_error
-- read_text
-- mkdir
-- file_to_write
-- list_dir
-- run_cmd
-- mp_wrapper
-- get_chunks
-- get_dups
-- concat_to_shm
-- concat_from_shm
-- revcomp
-- most_common
-- most_common_weighted
-- load_paths_txt
-- load_fasta
-- load_genbank
-
-Attributes:
------------
-- GZIP_EXT (str)
-- BASE_COMP (str.maketrans)
-- StartMethod (str, Enum)
 """
 
 __author__ = 'Michael X. Wang'
@@ -70,7 +37,7 @@
 BASE_COMP = str.maketrans('ATCGatcg', 'TAGCtagc') # Translation table for complement DNA bases
 
 class StartMethod(str, Enum):
-    """Start methods for multiprocessing. 
+    """Start methods for multiprocessing.
     """
     spawn = 'spawn'
     fork = 'fork'
@@ -83,12 +50,12 @@ class StartMethod(str, Enum):
 
 @dataclass(slots=True, frozen=True)
 class SharedArr:
-    """Class for a NumPy array attached to a SharedMemory instance. 
+    """Class for a NumPy array attached to a SharedMemory instance.
 
     Attributes:
-        name (str): Name of the SharedMemory instance. 
-        shape (tuple): Shape of the NumPy array. 
-        dtype (DTypeLike): dtype of the NumPy array. 
+        name (str): Name of the SharedMemory instance.
+        shape (tuple): Shape of the NumPy array.
+        dtype (DTypeLike): dtype of the NumPy array.
     """
     name: str
     shape: tuple[int, ...]
@@ -96,18 +63,18 @@ class SharedArr:
 
 
 def print_time_delta(seconds: float) -> None:
-    """Print time in seconds. 
+    """Print time in seconds.
     """
     logger.info(f' - Finished in {datetime.timedelta(seconds=seconds)}')
 
 
 def log_and_raise(
-    exception: type[Exception] = Exception, 
-    msg: str = '', 
-    from_none: bool = False, 
+    exception: type[Exception] = Exception,
+    msg: str = '',
+    from_none: bool = False,
     from_e: BaseException | None = None
 ) -> None:
-    """Log and raise an error. 
+    """Log and raise an error.
 
     Args:
         exception (type[Exception], optional): Exception type. [Exception]
@@ -129,29 +96,29 @@ def log_and_raise(
 
 
 def overwrite_warning(path: Path) -> None:
-    """Log overwrite warning. 
+    """Log overwrite warning.
     """
     logger.warning(f'File/directory already exists, content is overwritten (overwriting is turned on): {path}')
 
 
 def overwrite_error(path: Path) -> None:
-    """Raise FileExistsError. 
+    """Raise FileExistsError.
     """
     log_and_raise(FileExistsError, f'File/directory already exists, and overwriting is turned off: {path}', from_none=True)
 
 
 def read_text(path: Path) -> str:
-    """Read a UTF-8 text file with universal newline normalization. 
+    """Read a UTF-8 text file with universal newline normalization.
     """
     with open(path, 'r', encoding='utf-8', newline=None) as f:
         return f.read()
 
 
 def mkdir(path: Path, overwrite: bool=False, verbose=False) -> None:
-    """Create a directory. If the directory exists, remove it or raise an error. 
+    """Create a directory. If the directory exists, remove it or raise an error.
 
     Args:
-        path (Path): Path of the directory to be created. 
+        path (Path): Path of the directory to be created.
         overwrite (bool, optional): If True and the dir already exists, delete the existing dir and create a new empty one. [False]
         verbose (bool, optional): If True, print a warning message when content is overwritten. [False]
     """
@@ -170,10 +137,10 @@ def mkdir(path: Path, overwrite: bool=False, verbose=False) -> None:
 
 
 def file_to_write(path: Path, overwrite: bool=False, verbose=False) -> None:
-    """Prepare to write a file. If the file exists, remove it or raise an error. 
+    """Prepare to write a file. If the file exists, remove it or raise an error.
 
     Args:
-        path (Path): Path of the file to write. 
+        path (Path): Path of the file to write.
         overwrite (bool, optional): If True and the file already exists, delete the existing file. [False]
         verbose (bool, optional): If True, print a warning message when content is overwritten. [False]
     """
@@ -189,17 +156,17 @@ def file_to_write(path: Path, overwrite: bool=False, verbose=False) -> None:
 
 
 def list_dir(
-    path: Path = Path.cwd(), 
+    path: Path = Path.cwd(),
     mode: Literal['a', 'd', 'f'] = 'a'
 ) -> list[Path]:
-    """List all subdirectories and/or files under a directory. 
+    """List all subdirectories and/or files under a directory.
 
     Args:
         path (Path, optional): Path of the directory to list. ['./']
         mode (str, optional): 'd' to list subdirectories, 'f' to list files, 'a' to list all. ['a']
 
     Returns:
-        list: A list of sub-directories and/or files, sorted by names. 
+        list: A list of sub-directories and/or files, sorted by names.
     """
     # sanity check
     if not path.is_dir():
@@ -218,19 +185,19 @@ def list_dir(
 
 
 def run_cmd(
-    *args: str | Path, 
-    stdin: str | None = None, 
+    *args: str | Path,
+    stdin: str | None = None,
     raise_error: bool = True
 ) -> subprocess.CompletedProcess[str]:
-    """Run a command using subprocess.run(). Example usage: run_cmd('ls', '-a'). 
+    """Run a command using subprocess.run(). Example usage: run_cmd('ls', '-a').
 
     Args:
-        *args (str | Path): Command arguments. Must be strings or paths. 
+        *args (str | Path): Command arguments. Must be strings or paths.
         stdin (str, optional): Standard input. [None]
         raise_error (bool, optional): If True, raise an error if the command did not run successfully. [True]
 
     Returns:
-        CompletedProcess: command outputs, including stdout and stderr. 
+        CompletedProcess: command outputs, including stdout and stderr.
     """
     for a in args:
         if not isinstance(a, (str, Path)):
@@ -248,32 +215,32 @@ def run_cmd(
 
 
 def mp_wrapper(
-    func: Callable, 
-    all_args: Iterable, 
-    n_cpu: int=1, 
-    text: str | None=None, 
-    starmap: bool=True, 
-    unpack_output: bool=False, 
-    n_jobs: int | None=None, 
+    func: Callable,
+    all_args: Iterable,
+    n_cpu: int=1,
+    text: str | None=None,
+    starmap: bool=True,
+    unpack_output: bool=False,
+    n_jobs: int | None=None,
     start_method: StartMethod | None=_START_METHOD
 ) -> list:
-    """Wrapper for multiprocessing.Pool(). 
+    """Wrapper for multiprocessing.Pool().
 
     Args:
-        func (Callable): Function for multiprocessing. 
-        all_args (Iterable): Iterable of function arguments/parameters. 
+        func (Callable): Function for multiprocessing.
+        all_args (Iterable): Iterable of function arguments/parameters.
         n_cpu (int, optional): Number of processes to run in parallel. [1]
         text (str | None, optional): Message to be printed when multiprocessing starts. [None]
-        starmap (bool, optional): Use pool.starmap if True (func takes multiple arguments); 
+        starmap (bool, optional): Use pool.starmap if True (func takes multiple arguments);
             use pool.map if False (func takes only one argument). [True]
         unpack_output (bool, optional): If func has multiple output, return multiple lists instead of a single list of tuples. [False]
-        n_jobs (int | None, optional): Number of elements in `all_args`. 
+        n_jobs (int | None, optional): Number of elements in `all_args`.
             Helps determine the `chunksize` option for `pool.map` and `pool.starmap`. None to let Python decide. [None]
-        start_method (str | None, optional): Set the start method for multiprocessing ('fork', 'spawn', 'forkserver'). 
-            By default, 'spawn' is used for Windows and 'fork' for other systems. Use None to let Python decide. 
+        start_method (str | None, optional): Set the start method for multiprocessing ('fork', 'spawn', 'forkserver').
+            By default, 'spawn' is used for Windows and 'fork' for other systems. Use None to let Python decide.
 
     Returns:
-        list: A list of func outputs, in the same order as all_args. 
+        list: A list of func outputs, in the same order as all_args.
     """
     tik = time()
     if text:
@@ -312,14 +279,14 @@ def mp_wrapper(
 
 
 def get_chunks(ls: Sequence, n: int=1) -> Generator[Sequence, None, None]:
-    """Yield `n` roughly same size chunks of from a list. 
+    """Yield `n` roughly same size chunks of from a list.
 
     Args:
-        ls (Sequence): A list or list-alike. 
+        ls (Sequence): A list or list-alike.
         n (int): Number of chunks. [1]
 
     Yields:
-        Sequence: Chunks of ls. 
+        Sequence: Chunks of ls.
     """
     l = len(ls)
     size, remainder = divmod(l, n)
@@ -335,7 +302,7 @@ def get_chunks(ls: Sequence, n: int=1) -> Generator[Sequence, None, None]:
 
 
 def get_dups(iterable: Iterable[Hashable]) -> set:
-    """Returns a set of duplicated element(s) in an iterable. All elements should be Hashable. 
+    """Returns a set of duplicated element(s) in an iterable. All elements should be Hashable.
     """
     seen = set()
     duplicates = list()
@@ -348,15 +315,15 @@ def get_dups(iterable: Iterable[Hashable]) -> set:
 
 
 def concat_to_shm(arrays: Sequence[NDArray]) -> SharedArr:
-    """Concat NumPy arrays along the first dimension into a shared memory block. 
-    - Each individual array is deleted during this process to save memory. 
-    - Arrays are copied as raw bytes to bypass overhead of NumPy assignment. 
+    """Concat NumPy arrays along the first dimension into a shared memory block.
+    - Each individual array is deleted during this process to save memory.
+    - Arrays are copied as raw bytes to bypass overhead of NumPy assignment.
 
     Args:
-        arrays (Sequence[NDArray]): Arrays to be concatenated. 
+        arrays (Sequence[NDArray]): Arrays to be concatenated.
 
     Returns:
-        SharedArr: The concatenated array attached to a SharedMemory instance. 
+        SharedArr: The concatenated array attached to a SharedMemory instance.
     """
     if not arrays:
         log_and_raise(ValueError, 'No array is provided')
@@ -376,7 +343,7 @@ def concat_to_shm(arrays: Sequence[NDArray]) -> SharedArr:
     n0 = sum(arr.shape[0] for arr in arrays)
     out_shape = (n0, *trailing_shape)
     out_shm = SharedMemory(
-        create=True, 
+        create=True,
         size=int(np.prod(out_shape, dtype=np.int64) * dtype.itemsize)
     )
     try:
@@ -402,16 +369,16 @@ def concat_to_shm(arrays: Sequence[NDArray]) -> SharedArr:
 
 
 def concat_from_shm(arrays: Sequence[SharedArr], n_cpu: int=1) -> NDArray:
-    """Concat SharedArr instances along the first dimension into a NumPy array. 
-    - Each SharedArr is unlinked during this process to save memory. 
-    - Arrays are copied as raw bytes, in parallel. 
+    """Concat SharedArr instances along the first dimension into a NumPy array.
+    - Each SharedArr is unlinked during this process to save memory.
+    - Arrays are copied as raw bytes, in parallel.
 
     Args:
-        arrays (Sequence[SharedArr]): Arrays to be concatenated. 
+        arrays (Sequence[SharedArr]): Arrays to be concatenated.
         n_cpu (int, optional): Number of threads to run in parallel. [1]
 
     Returns:
-        NDArray: The concatenated NumPy array. 
+        NDArray: The concatenated NumPy array.
     """
     if not arrays:
         log_and_raise(ValueError, 'No array is provided')
@@ -437,7 +404,7 @@ def concat_from_shm(arrays: Sequence[SharedArr], n_cpu: int=1) -> NDArray:
     # copy each array into its slot in out_buf as raw bytes, in parallel
     # we have to use thread-based parallelism to write to the same private memory block
     def copy_array(arr_name: str, offset: int, arr_size: int):
-        """The worker function for each thread. 
+        """The worker function for each thread.
         """
         arr_shm = SharedMemory(arr_name)
         try:
@@ -469,35 +436,35 @@ def copy_array(arr_name: str, offset: int, arr_size: int):
 
 
 def revcomp(seq: str) -> str:
-    """Return the reverse complement of a DNA sequence. 
+    """Return the reverse complement of a DNA sequence.
     """
     return seq.translate(BASE_COMP)[::-1]
 
 
 def most_common(iterable: Iterable[Hashable]):
-    """Return the most common element in an Iterable. 
-    All elements should be Hashable. 
+    """Return the most common element in an Iterable.
+    All elements should be Hashable.
     """
     return Counter(iterable).most_common(1)[0][0]
 
 
 def most_common_weighted(iterable: Iterable):
-    """Return the most common element in an Iterable, weighted by element length. 
-    Each element should be Hashable and Sized (e.g., tuple or str). 
+    """Return the most common element in an Iterable, weighted by element length.
+    Each element should be Hashable and Sized (e.g., tuple or str).
     """
     c = Counter(iterable)
     return max(c, key=lambda k: len(k)*c[k])
 
 
 def load_paths_txt(paths_txt: Path) -> list[Path]:
-    """Load file paths from a text file, with one path per line. 
-    Relative paths are resolved relative to the directory containing `paths_txt`. 
+    """Load file paths from a text file, with one path per line.
+    Relative paths are resolved relative to the directory containing `paths_txt`.
 
     Args:
-        paths_txt (Path): A text file with one path per line. 
+        paths_txt (Path): A text file with one path per line.
 
     Returns:
-        list[Path]: A list of valid and resolved file paths. 
+        list[Path]: A list of valid and resolved file paths.
     """
     paths_txt = paths_txt.resolve(strict=True)
     base_dir = paths_txt.parent
@@ -523,14 +490,14 @@ def load_paths_txt(paths_txt: Path) -> list[Path]:
 
 
 def load_fasta(path: Path) -> tuple[str, ...]:
-    """Parse an assembly file in FASTA format and return its sequences. 
-    Gzip files are supported (file name should end with .gz). 
+    """Parse an assembly file in FASTA format and return its sequences.
+    Gzip files are supported (file name should end with .gz).
 
     Args:
-        path (Path): Path to the FASTA file. If the file is gzipped, the extension should be .gz. 
+        path (Path): Path to the FASTA file. If the file is gzipped, the extension should be .gz.
 
     Returns:
-        tuple[str, ...]: Sequences of FASTA records (upper case), in the same order as they appear in the file. 
+        tuple[str, ...]: Sequences of FASTA records (upper case), in the same order as they appear in the file.
     """
     # read file content
     if path.suffix == GZIP_EXT:
@@ -565,16 +532,16 @@ def load_fasta(path: Path) -> tuple[str, ...]:
 
 if _HAS_BIO:
     def load_genbank(path: Path) -> tuple[dict[str, str], int]:
-        """Parse an assembly file in GenBank format, and get the sequences and IDs of all records. 
-        Gzip files are supported (file name should end with .gz). 
+        """Parse an assembly file in GenBank format, and get the sequences and IDs of all records.
+        Gzip files are supported (file name should end with .gz).
 
         Args:
-            path (Path): Path to the GenBank file. If the file is gzipped, the extension should be .gz. 
+            path (Path): Path to the GenBank file. If the file is gzipped, the extension should be .gz.
 
         Returns:
             tuple: A tuple containing
-                1. dict[str, str]: A dictionary with record ID as key and record sequence as value. 
-                2. int: Sum of the length of all sequence records. 
+                1. dict[str, str]: A dictionary with record ID as key and record sequence as value.
+                2. int: Sum of the length of all sequence records.
         """
         if path.suffix == GZIP_EXT:
             with gzip.open(path, 'rb') as f:
@@ -592,7 +559,7 @@ def load_genbank(path: Path) -> tuple[dict[str, str], int]:
             all_record[record_id] = seq
             all_id.append(record_id)
             total_len += len(seq)
-        
+
         # check duplicate record ID
         if len(all_id) != len(set(all_id)):
             logger.warning(f' - Duplicate record ID(s) {get_dups(all_id)}, in: {path}')
@@ -600,7 +567,7 @@ def load_genbank(path: Path) -> tuple[dict[str, str], int]:
 else:
     def load_genbank(path) -> None:
         log_and_raise(
-            ImportError, 
-            'Biopython is needed for parsing GenBank files', 
+            ImportError,
+            'Biopython is needed for parsing GenBank files',
             from_none=True
         )
\ No newline at end of file
diff --git a/tests/smoke/test_cli.py b/tests/smoke/test_cli.py
index 83becb0..9bc7b00 100644
--- a/tests/smoke/test_cli.py
+++ b/tests/smoke/test_cli.py
@@ -1,16 +1,21 @@
+import re
 from pathlib import Path
 
-from click import unstyle
 from typer.testing import CliRunner
 
 from seqwin import __version__
 from seqwin import cli
 
+ANSI_RE = re.compile(r"\x1b\[[0-?]*[ -/]*[@-~]")
 runner = CliRunner()
 
 
+def unstyle(text: str) -> str:
+    return ANSI_RE.sub("", text)
+
+
 def test_help_shows_key_options() -> None:
-    result = runner.invoke(cli.app, ['--help'], terminal_width=120)
+    result = runner.invoke(cli.app, ['--help'], terminal_width=120, color=False)
     output = unstyle(result.output)
 
     assert result.exit_code == 0
diff --git a/tests/smoke/test_graph.py b/tests/smoke/test_graph.py
index 18af0a5..f77e076 100644
--- a/tests/smoke/test_graph.py
+++ b/tests/smoke/test_graph.py
@@ -48,31 +48,31 @@ def test_build_threading_equivalence(targets_dir, non_targets_dir) -> None:
         non_targets_dir / 'non-target-1.fasta',
         non_targets_dir / 'non-target-2.fasta',
     ]
-    assembly_idx = [0, 1, 2, 3]
-    is_target = [True, True, False, False]
+    assembly_indices = [0, 1, 2, 3]
+    is_targets = [True, True, False, False]
 
     kmers_1, idx_1, nodes_1, edges_1, record_ids_1 = build(
         assembly_paths,
         kmerlen=7,
         windowsize=10,
-        assembly_idx=assembly_idx,
-        is_target=is_target,
+        assembly_indices=assembly_indices,
+        is_targets=is_targets,
         n_cpu=1,
     )
     kmers_2, idx_2, nodes_2, edges_2, record_ids_2 = build(
         assembly_paths,
         kmerlen=7,
         windowsize=10,
-        assembly_idx=assembly_idx,
-        is_target=is_target,
+        assembly_indices=assembly_indices,
+        is_targets=is_targets,
         n_cpu=2,
     )
     kmers_many, idx_many, nodes_many, edges_many, record_ids_many = build(
         assembly_paths,
         kmerlen=7,
         windowsize=10,
-        assembly_idx=assembly_idx,
-        is_target=is_target,
+        assembly_indices=assembly_indices,
+        is_targets=is_targets,
         n_cpu=99,
     )