spoa_rs 0.1.0

Rust bindings for SPOA (SIMD POA)
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

// Copyright (c) 2020 Robert Vaser

#include "spoa/alignment_engine.hpp"

#include <algorithm>
#include <exception>
#include <limits>
#include <stdexcept>

#include "sisd_alignment_engine.hpp"
#include "simd_alignment_engine.hpp"

namespace spoa {

std::unique_ptr<AlignmentEngine> AlignmentEngine::Create(
    AlignmentType type,
    std::int8_t m,
    std::int8_t n,
    std::int8_t g) {
  return Create(type, m, n, g, g);
}

std::unique_ptr<AlignmentEngine> AlignmentEngine::Create(
    AlignmentType type,
    std::int8_t m,
    std::int8_t n,
    std::int8_t g,
    std::int8_t e) {
  return Create(type, m, n, g, e, g, e);
}

std::unique_ptr<AlignmentEngine> AlignmentEngine::Create(
    AlignmentType type,
    std::int8_t m,
    std::int8_t n,
    std::int8_t g,
    std::int8_t e,
    std::int8_t q,
    std::int8_t c) {
  if (type != AlignmentType::kSW &&
      type != AlignmentType::kNW &&
      type != AlignmentType::kOV) {
    throw std::invalid_argument(
        "[spoa::AlignmentEngine::Create] error: invalid alignment type!");
  }
  if (g > 0 || q > 0) {
    throw std::invalid_argument(
        "[spoa::AlignmentEngine::Create] error: "
        "gap opening penalty must be non-positive!");
  }
  if (e > 0 || c > 0) {
    throw std::invalid_argument(
        "[spoa::AlignmentEngine::Create] error: "
        "gap extension penalty must be non-positive!");
  }

  AlignmentSubtype subtype = g >= e ?
      AlignmentSubtype::kLinear : (g <= q || e >= c ?
      AlignmentSubtype::kAffine : AlignmentSubtype::kConvex);

  if (subtype == AlignmentSubtype::kLinear) {
    e = g;
  } else if (subtype == AlignmentSubtype::kAffine) {
    q = g;
    c = e;
  }

  auto dst = CreateSimdAlignmentEngine(type, subtype, m, n, g, e, q, c);
  if (!dst) {
    return SisdAlignmentEngine::Create(type, subtype, m, n, g, e, q, c);
  }
  return dst;
}

AlignmentEngine::AlignmentEngine(
    AlignmentType type,
    AlignmentSubtype subtype,
    std::int8_t m,
    std::int8_t n,
    std::int8_t g,
    std::int8_t e,
    std::int8_t q,
    std::int8_t c)
    : type_(type),
      subtype_(subtype),
      m_(m),
      n_(n),
      g_(g),
      e_(e),
      q_(q),
      c_(c) {
}

Alignment AlignmentEngine::Align(
    const std::string& sequence,
    const Graph& graph,
    std::int32_t* score) {
  return Align(sequence.c_str(), sequence.size(), graph, score);
}

std::int64_t AlignmentEngine::WorstCaseAlignmentScore(
    std::int64_t i,
    std::int64_t j) const {
  auto gap_score = [&] (std::int64_t len) -> std::int64_t {
    return len == 0 ? 0 : std::min(g_ + (len - 1) * e_, q_ + (len - 1) * c_);
  };
  return std::min(
      -1 * (m_ * std::min(i, j) + gap_score(std::abs(i - j))),
      gap_score(i) + gap_score(j));
}

}  // namespace spoa