Updates to replace pointers with iterators for BAM and SAM inputs

CMakeLists.txt

@@ -104,6 +104,7 @@ target_link_libraries(fastq_bgzf_file PUBLIC
 
 add_library(bam_header OBJECT src/bam_header.cpp)
 
+add_library(bamrec OBJECT src/bamrec.cpp)
 add_library(bam_file OBJECT src/bam_file.cpp)
 target_link_libraries(bam_file PUBLIC HTSLIB::HTSLIB)
 
@@ -210,6 +211,7 @@ target_link_libraries(falco PUBLIC
   fastq_gz_file
   fastq_bgzf_file
   bam_file
+  bamrec
   sam_file
   samrec
   bam_header

src/bam_file.cpp

@@ -7,6 +7,7 @@
 #include "falco_word.hpp"
 #include "task_queue.hpp"
 
+#include <htslib/bgzf.h>  // IWYU pragma: keep
 #include <htslib/hfile.h>
 #include <htslib/sam.h>
 
@@ -40,7 +41,7 @@ estimate_n_reads_bam(const std::string &filename)
     throw std::runtime_error("failed to identify file format: " + filename);
 
   // NOLINTNEXTLINE(cppcoreguidelines-pro-type-union-access)
-  const auto &fp = format->format == bam ? f->fp.bgzf->fp : f->fp.hfile;
+  const auto fp = format->format == bam ? f->fp.bgzf->fp : f->fp.hfile;
   const auto pos_after_header = htell(fp);
   std::unique_ptr<bam1_t, void (*)(bam1_t *)> rec(bam_init1(), &bam_destroy1);
   std::uint64_t n_reads{};
@@ -53,11 +54,12 @@ estimate_n_reads_bam(const std::string &filename)
     throw std::system_error(std::make_error_code(std::errc(errno)),
                             "error reading bam record from: " + filename);
   const auto pos_after_reads = htell(fp);
-  const auto n_compressed_bytes = pos_after_reads - pos_after_header;
+  const auto n_bytes = pos_after_reads - pos_after_header;
   const auto filesize = std::filesystem::file_size(filename);
-  const auto estimate = static_cast<std::uint64_t>(
-    as_frac(n_reads * (filesize - pos_after_header), n_compressed_bytes));
-  return {estimate, total_read_len / n_reads, filesize};
+  const auto n_reads_estimate = static_cast<std::uint64_t>(
+    as_frac(n_reads * (filesize - pos_after_header), n_bytes));
+  const auto read_len_estimate = total_read_len / n_reads;
+  return {n_reads_estimate, read_len_estimate, filesize};
 }
 
 [[nodiscard]] auto
@@ -138,10 +140,10 @@ partition(auto itr,                     //
           std::atomic_int32_t &n_tasks) {
   // ADS: this isn't working as desired: the end position of each part should be
   // the first record end past the 'end_itr' below unless end_itr == end
-  const auto dist = std::distance(itr, end);
+  auto dist = std::distance(itr, end);
   const auto chunk_size = (dist + n_chunks - 1) / n_chunks;
   while (itr != end) {
-    const auto dist = std::distance(itr, end);
+    dist = std::distance(itr, end);
     auto end_itr = itr + (dist < chunk_size ? dist : chunk_size);
     // ADS: find_end_pos doesn't find end pos of a record, but of a range, so
     // includes multiple records
@@ -164,13 +166,14 @@ bam_file::get_chunks(const std::int64_t n_chunks,  //
   // input buffer has been inflated and will provide data for analysis
   std::swap(input_buffer, output_buffer);
   std::swap(output_last, input_last);
-  auto itr = std::data(output_buffer);
-  const auto end = itr + output_last;  // NOLINT(*-pointer-arithmetic)
-  if (bh)                              // if we are still parsing the header
+  const auto beg = std::cbegin(output_buffer);
+  const auto end = beg + output_last;
+  auto itr = beg;
+  if (bh)  // if we are still parsing the header
     itr = bh.update(itr, end);
-  if (std::distance(itr, end) > 0)
+  if (itr != end)
     itr = partition(itr, end, n_chunks, file_id, tq, n_tasks);
-  output_cursor = std::distance(std::data(output_buffer), itr);
+  output_cursor = std::distance(beg, itr);
 }
 
 auto

src/bam_header.cpp

@@ -6,8 +6,8 @@
 #include <stdexcept>
 
 [[nodiscard]] auto
-bam_header::update(bam_header::iterator itr,
-                   const bam_header::iterator end) -> bam_header::iterator {
+bam_header::update(const_iterator itr,
+                   const const_iterator end) -> const_iterator {
   static constexpr auto msg = "incorrect BAM magic identified: {} at {}";
   const auto update_u32 = [](auto &val, const auto inc, const auto the_byte) {
     // NOLINTNEXTLINE(*-avoid-magic-numbers)

src/bam_header.hpp

@@ -7,10 +7,11 @@
 
 #include <cstdint>
 #include <string>
+#include <vector>
 
 struct bam_header {
-  using iterator = char *;
-  using const_iterator = const char *;
+  using iterator = std::vector<char>::iterator;
+  using const_iterator = std::vector<char>::const_iterator;
 
   static constexpr auto magic = "BAM\1";
   static constexpr auto magic_size = 4;
@@ -30,7 +31,7 @@ struct bam_header {
   }
 
   [[nodiscard]] auto
-  update(iterator itr, const iterator end) -> iterator;
+  update(const_iterator itr, const const_iterator end) -> const_iterator;
 
   [[nodiscard]] auto
   ref_incomplete() const -> bool {

src/bamrec.cpp

@@ -0,0 +1,121 @@
+// SPDX-License-Identifier: MIT; Copyright 2026 Andrew D Smith
+
+#include "bamrec.hpp"
+
+#include <algorithm>
+#include <cstring>
+#include <format>
+#include <memory>
+#include <span>
+#include <string>
+#include <utility>
+
+[[nodiscard]] auto
+bamrec::to_string() const -> std::string {
+  const auto buffer_s = std::span(std::cbegin(buffer), std::cend(buffer));
+  const auto name_itr = std::cbegin(buffer_s);
+  const auto seq_itr = name_itr + name_len;
+  const auto qual_itr = seq_itr + seq_len;
+  std::string qual_fixed(seq_len, '\0');
+  std::transform(qual_itr, qual_itr + seq_len, std::begin(qual_fixed),
+                 [](const auto c) { return c + quality_score_offset; });
+  return std::format("@{}\n{}\n+\n{}",                //
+                     std::string(name_itr, seq_itr),  //
+                     std::string(seq_itr, qual_itr),  //
+                     qual_fixed);
+}
+
+#ifdef seq_nt16_str
+#undef seq_nt16_str
+#endif
+
+#ifdef bam_seqi
+#undef bam_seqi
+#endif
+
+template <class BidirIt, class OutputIt>
+static inline constexpr OutputIt
+assign_sequence_revcomp(BidirIt first, auto last, OutputIt d_first) {
+  constexpr auto complem = [](const auto x) {
+    return "TNGNNNCNNNNNNNNNNNNA"[x - 'A'];
+  };
+  constexpr auto seq_nt16_str = "=ACMGRSVTWYHKDBN";
+  constexpr auto bam_seqi = [](const auto s, const auto i) -> int {
+    constexpr auto low_nibble_on = 0xf;
+    return s[i >> 1] >> ((~i & 1) << 2) & low_nibble_on;
+  };
+  for (auto j = last; j != 0; ++d_first)
+    // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
+    *d_first = complem(seq_nt16_str[bam_seqi(first, --j)]);
+  return d_first;
+}
+
+template <class BidirIt, class OutputIt>
+static inline constexpr OutputIt
+assign_sequence(BidirIt first, auto last, OutputIt d_first) {
+  constexpr auto seq_nt16_str = "=ACMGRSVTWYHKDBN";
+  constexpr auto bam_seqi = [](const auto s, const auto i) -> int {
+    constexpr auto low_nibble_on = 0xf;
+    return s[i >> 1] >> ((~i & 1) << 2) & low_nibble_on;
+  };
+  for (auto j = 0U; j != last; ++j)
+    // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
+    *d_first++ = seq_nt16_str[bam_seqi(first, j)];
+  return d_first;
+}
+
+[[nodiscard]] auto
+bamrec::get_next(bamrec::pos_t &itr, const bamrec::pos_t end,
+                 bamrec &rec) -> bool {
+  if (std::distance(itr, end) < bam_core_t::sz)
+    return false;
+  bam_core_t core{};
+  std::memcpy(std::addressof(core), std::to_address(itr), sizeof(bam_core_t));
+  if (std::distance(itr, end) < core.real_block_size())
+    return false;
+  rec.name_len = core.l_read_name - 1;  // we don't need the '\0'
+  rec.seq_len = core.l_seq;
+  const auto rec_size = rec.name_len + 2 * rec.seq_len;
+  if (std::size(rec.buffer) < rec_size)
+    rec.buffer.resize(rec_size);
+  auto out_itr = std::begin(rec.buffer);
+  std::copy_n(itr + bam_core_t::read_name_offset, rec.name_len, out_itr);
+  out_itr += rec.name_len;  // increment data cursor to sequence
+  auto seq_in = itr + core.seq_offset();
+  if (core.bam_is_rev())
+    assign_sequence_revcomp(seq_in, core.l_seq, out_itr);
+  else
+    assign_sequence(seq_in, core.l_seq, out_itr);
+  out_itr += rec.seq_len;  // increment data cursor to qual
+  const auto has_qual = (itr[core.qual_offset()] != qual_missing_code);
+  if (has_qual) {
+    const auto qual_itr = itr + core.qual_offset();
+    const auto qual_end = qual_itr + rec.seq_len;
+    if (core.bam_is_rev())
+      std::reverse_copy(qual_itr, qual_end, out_itr);
+    else
+      std::copy(qual_itr, qual_end, out_itr);
+  }
+  else
+    *out_itr = static_cast<char>(qual_missing_code);
+  itr += core.real_block_size();
+  return true;
+}
+
+[[nodiscard]] auto
+bamrec::find_end_pos(bamrec::pos_t itr,
+                     const bamrec::pos_t end) -> bamrec::pos_t {
+  static constexpr std::int64_t record_size_size = sizeof(std::uint32_t);
+  std::uint32_t record_size{};
+  while (itr != end) {
+    if (std::distance(itr, end) < record_size_size)
+      return itr;
+    std::memcpy(std::addressof(record_size), std::to_address(itr),
+                record_size_size);
+    record_size += record_size_size;
+    if (std::distance(itr, end) < record_size)
+      return itr;
+    itr += record_size;  // only increment on consume of full record
+  }
+  return itr;
+}