bihash_kmer_directional.hpp

#pragma once
 
#include <array>
#include <cstdint>
#include <cstring>
#include <deque>
#include <memory>
#include <string>
#include <vector>
 
#include <btllib/bihash_hashing_internals.hpp>
#include <btllib/hashing_internals.hpp>
#include <btllib/nthash_kmer.hpp>
#include <btllib/status.hpp>
 
namespace btllib {
 
using hashing_internals::base_forward_hash;
using hashing_internals::base_reverse_hash;
using hashing_internals::extend_hashes;
using hashing_internals::next_forward_hash;
using hashing_internals::next_reverse_hash;
using hashing_internals::prev_forward_hash;
using hashing_internals::prev_reverse_hash;
using hashing_internals::SEED_N;
using hashing_internals::SEED_TAB;
using hashing_internals::sub_hash;
 
using hashing_internals::BI_CT_SEED_TAB;
using hashing_internals::BI_GA_SEED_TAB;
 
using hashing_internals::BI_CT_MS_TAB_33R;
using hashing_internals::BI_GA_MS_TAB_33R;
 
using hashing_internals::BI_CT_MS_TAB_31L;
using hashing_internals::BI_GA_MS_TAB_31L;
 
using hashing_internals::BI_CT_CONVERT_TAB;
using hashing_internals::BI_GA_CONVERT_TAB;
 
using hashing_internals::BI_CT_RC_CONVERT_TAB;
using hashing_internals::BI_GA_RC_CONVERT_TAB;
 
using hashing_internals::BI_CT_DIMER_TAB;
using hashing_internals::BI_CT_TETRAMER_TAB;
using hashing_internals::BI_CT_TRIMER_TAB;
using hashing_internals::BI_GA_DIMER_TAB;
using hashing_internals::BI_GA_TETRAMER_TAB;
using hashing_internals::BI_GA_TRIMER_TAB;
 
class BiHashDirectional
{
 
public:
  BiHashDirectional(const char* seq,
                    size_t seq_len,
                    hashing_internals::NUM_HASHES_TYPE num_hashes,
                    hashing_internals::K_TYPE k,
                    std::string conversion_type,
                    size_t pos = 0)
    : seq(seq)
    , seq_len(seq_len)
    , num_hashes(num_hashes)
    , k(k)
    , conversion_type(conversion_type)
    , pos(pos)
    , initialized(false)
    , hash_arr(new uint64_t[num_hashes])
  {
    check_error(k == 0, "BiHash: k must be greater than 0");
    check_error(k % 2 == 1, "BiHash: k must be even");
    check_error((k / 2) % 2 == 0, "BiHash: k must have odd number of dimers");
    check_error(this->seq_len < k,
                "BiHash: sequence length (" + std::to_string(this->seq_len) +
                  ") is smaller than k (" + std::to_string(k) + ")");
    check_error(pos > this->seq_len - k,
                "BiHash: passed position (" + std::to_string(pos) +
                  ") is larger than sequence length (" +
                  std::to_string(this->seq_len) + ")");
 
    size_t num_dimers = k / 2;
    size_t center_dimer_start = (num_dimers / 2) * 2;
    if (conversion_type == "CT") {
      meth_base_idx = center_dimer_start;
 
      primitive_tab = BI_CT_SEED_TAB;
 
      right_table = BI_CT_MS_TAB_33R;
      left_table = BI_CT_MS_TAB_31L;
 
      convert_tab = BI_CT_CONVERT_TAB;
      rc_convert_tab = BI_CT_RC_CONVERT_TAB;
 
      dimer_tab = BI_CT_DIMER_TAB;
      rc_dimer_tab = BI_GA_DIMER_TAB;
      trimer_tab = BI_CT_TRIMER_TAB;
      rc_trimer_tab = BI_GA_TRIMER_TAB;
      tetramer_tab = BI_CT_TETRAMER_TAB;
      ;
      rc_tetramer_tab = BI_GA_TETRAMER_TAB;
    } else if (conversion_type == "GA") {
      meth_base_idx = center_dimer_start + 1;
 
      primitive_tab = BI_GA_SEED_TAB;
 
      right_table = BI_GA_MS_TAB_33R;
      left_table = BI_GA_MS_TAB_31L;
 
      convert_tab = BI_GA_CONVERT_TAB;
      rc_convert_tab = BI_GA_RC_CONVERT_TAB;
 
      dimer_tab = BI_GA_DIMER_TAB;
      rc_dimer_tab = BI_CT_DIMER_TAB;
      trimer_tab = BI_GA_TRIMER_TAB;
      rc_trimer_tab = BI_CT_TRIMER_TAB;
      tetramer_tab = BI_GA_TETRAMER_TAB;
      rc_tetramer_tab = BI_CT_TETRAMER_TAB;
    } else {
      check_error(
        true, "Unsupported conversion_type. Must be CT or GA, as capitalized.");
    }
  }
 
  BiHashDirectional(const std::string& seq,
                    hashing_internals::NUM_HASHES_TYPE num_hashes,
                    hashing_internals::K_TYPE k,
                    std::string conversion_type,
                    size_t pos = 0)
    : BiHashDirectional(seq.data(),
                        seq.size(),
                        num_hashes,
                        k,
                        conversion_type,
                        pos)
  {
  }
 
  BiHashDirectional(const BiHashDirectional& obj)
    : seq(obj.seq)
    , seq_len(obj.seq_len)
    , num_hashes(obj.num_hashes)
    , k(obj.k)
    , conversion_type(obj.conversion_type)
    , pos(obj.pos)
    , initialized(obj.initialized)
    , fwd_hash(obj.fwd_hash)
    , rev_hash(obj.rev_hash)
    , hash_arr(new uint64_t[obj.num_hashes])
    , primitive_tab(obj.primitive_tab)
    , right_table(obj.right_table)
    , left_table(obj.left_table)
    , convert_tab(obj.convert_tab)
    , rc_convert_tab(obj.rc_convert_tab)
    , dimer_tab(obj.dimer_tab)
    , rc_dimer_tab(obj.rc_dimer_tab)
    , trimer_tab(obj.trimer_tab)
    , rc_trimer_tab(obj.rc_trimer_tab)
    , tetramer_tab(obj.tetramer_tab)
    , rc_tetramer_tab(obj.rc_tetramer_tab)
  {
    std::memcpy(
      hash_arr.get(), obj.hash_arr.get(), num_hashes * sizeof(uint64_t));
  }
 
  BiHashDirectional(BiHashDirectional&&) = default;
  bool roll()
  {
    if (!initialized) {
      return init();
    }
    if (pos >= seq_len - k) {
      return false;
    }
    if (hashing_internals::SEED_TAB[(unsigned char)seq[pos + k]] ==
        hashing_internals::SEED_N) {
      pos += k;
      return init();
    }
    fwd_hash = next_forward_hash(fwd_hash,
                                 k,
                                 seq[pos],
                                 seq[pos + k],
                                 primitive_tab,
                                 right_table,
                                 left_table);
    rev_hash = next_reverse_hash(rev_hash,
                                 k,
                                 seq[pos],
                                 seq[pos + k],
                                 primitive_tab,
                                 right_table,
                                 left_table);
    extend_hashes(fwd_hash, rev_hash, k, num_hashes, hash_arr.get());
    ++pos;
    return true;
  }
 
  bool roll_back()
  {
    if (!initialized) {
      return init();
    }
    if (pos == 0) {
      return false;
    }
    if (SEED_TAB[(unsigned char)seq[pos - 1]] == SEED_N && pos >= k) {
      pos -= k;
      return init();
    }
    if (SEED_TAB[(unsigned char)seq[pos - 1]] == SEED_N) {
      return false;
    }
    fwd_hash = prev_forward_hash(fwd_hash,
                                 k,
                                 seq[pos + k - 1],
                                 seq[pos - 1],
                                 primitive_tab,
                                 right_table,
                                 left_table);
    rev_hash = prev_reverse_hash(rev_hash,
                                 k,
                                 seq[pos + k - 1],
                                 seq[pos - 1],
                                 primitive_tab,
                                 right_table,
                                 left_table);
    extend_hashes(fwd_hash, rev_hash, k, num_hashes, hash_arr.get());
    --pos;
    return true;
  }
 
  bool peek()
  {
    if (pos >= seq_len - k) {
      return false;
    }
    return peek(seq[pos + k]);
  }
 
  bool peek_back()
  {
    if (pos == 0) {
      return false;
    }
    return peek_back(seq[pos - 1]);
  }
 
  bool peek(char char_in)
  {
    if (!initialized) {
      return init();
    }
    if (SEED_TAB[(unsigned char)char_in] == SEED_N) {
      return false;
    }
    const uint64_t fwd = next_forward_hash(
      fwd_hash, k, seq[pos], char_in, primitive_tab, right_table, left_table);
    const uint64_t rev = next_reverse_hash(
      rev_hash, k, seq[pos], char_in, primitive_tab, right_table, left_table);
    extend_hashes(fwd, rev, k, num_hashes, hash_arr.get());
    return true;
  }
 
  bool peek_back(char char_in)
  {
    if (!initialized) {
      return init();
    }
    if (SEED_TAB[(unsigned char)char_in] == SEED_N) {
      return false;
    }
    const unsigned char char_out = seq[pos + k - 1];
    const uint64_t fwd = prev_forward_hash(
      fwd_hash, k, char_out, char_in, primitive_tab, right_table, left_table);
    const uint64_t rev = prev_reverse_hash(
      rev_hash, k, char_out, char_in, primitive_tab, right_table, left_table);
    extend_hashes(fwd, rev, k, num_hashes, hash_arr.get());
    return true;
  }
 
  bool is_methylated()
  {
    char base = seq[pos + meth_base_idx];
 
    if (conversion_type == "CT")
      return base == 'C' || base == 'c';
    else
      return base == 'G' || base == 'g';
  }
 
  bool is_forward_higher() const { return fwd_hash > rev_hash; }
 
  void sub(const std::vector<unsigned>& positions,
           const std::vector<unsigned char>& new_bases)
  {
    sub_hash(fwd_hash,
             rev_hash,
             seq + pos,
             positions,
             new_bases,
             get_k(),
             get_hash_num(),
             hash_arr.get(),
             right_table,
             left_table);
  }
 
  const uint64_t* hashes() const { return hash_arr.get(); }
 
  size_t get_pos() const { return pos; }
 
  hashing_internals::NUM_HASHES_TYPE get_hash_num() const { return num_hashes; }
 
  hashing_internals::K_TYPE get_k() const { return k; }
 
  uint64_t get_forward_hash() const { return fwd_hash; }
 
  uint64_t get_reverse_hash() const { return rev_hash; }
 
private:
  const char* seq;
  const size_t seq_len;
  hashing_internals::NUM_HASHES_TYPE num_hashes;
  hashing_internals::K_TYPE k;
  std::string conversion_type;
  size_t pos;
  size_t meth_base_idx;
  bool initialized;
  uint64_t fwd_hash = 0;
  uint64_t rev_hash = 0;
  std::unique_ptr<uint64_t[]> hash_arr;
 
  const uint64_t* primitive_tab;
 
  const uint64_t* const* right_table;
  const uint64_t* const* left_table;
 
  const uint8_t* convert_tab;
  const uint8_t* rc_convert_tab;
 
  const uint64_t* dimer_tab;
  const uint64_t* rc_dimer_tab;
  const uint64_t* trimer_tab;
  const uint64_t* rc_trimer_tab;
  const uint64_t* tetramer_tab;
  const uint64_t* rc_tetramer_tab;
 
  bool init()
  {
    bool has_n = true;
    while (pos <= seq_len - k + 1 && has_n) {
      has_n = false;
      for (unsigned i = 0; i < k && pos <= seq_len - k + 1; i++) {
        if (SEED_TAB[(unsigned char)seq[pos + k - i - 1]] == SEED_N) {
          pos += k - i;
          has_n = true;
        }
      }
    }
    if (pos > seq_len - k) {
      return false;
    }
    fwd_hash = base_forward_hash(seq + pos,
                                 k,
                                 primitive_tab,
                                 dimer_tab,
                                 trimer_tab,
                                 tetramer_tab,
                                 convert_tab);
    rev_hash = base_reverse_hash(seq + pos,
                                 k,
                                 primitive_tab,
                                 rc_dimer_tab,
                                 rc_trimer_tab,
                                 rc_tetramer_tab,
                                 rc_convert_tab);
    extend_hashes(fwd_hash, rev_hash, k, num_hashes, hash_arr.get());
    initialized = true;
    return true;
  }
};
 
} // namespace btllib