ScanTemplate.hpp

Go to the documentation of this file.

#ifndef KAORI_SCAN_TEMPLATE_HPP
#define KAORI_SCAN_TEMPLATE_HPP
 
#include <bitset>
#include <vector>
#include <stdexcept>
#include <string>
 
#include "utils.hpp"
 
namespace kaori {
 
template<SeqLength max_size_>
class ScanTemplate { 
private:
    static constexpr SeqLength N = max_size_ * 4;
 
public:
    ScanTemplate() = default;
 
    ScanTemplate(const char* template_seq, SeqLength template_length, SearchStrand strand) :
        my_length(template_length),
        my_forward(search_forward(strand)),
        my_reverse(search_reverse(strand))
    {
        if (my_length > max_size_) {
            throw std::runtime_error("maximum template size should be " + std::to_string(max_size_) + " bp");
        }
 
        if (my_forward) {
            for (SeqLength i = 0; i < my_length; ++i) {
                char b = template_seq[i];
                if (b != '-') {
                    add_base_to_hash(my_forward_ref, b);
                    add_mask_to_hash(my_forward_mask);
                    my_forward_mask_ambiguous.set(i);
                } else {
                    shift_hash(my_forward_ref);
                    shift_hash(my_forward_mask);
                    add_variable_base(my_forward_variables, i);
                }
            }
        } else {
            // Forward variable regions are always defined.
            for (SeqLength i = 0; i < my_length; ++i) {
                char b = template_seq[i];
                if (b == '-') {
                    add_variable_base(my_forward_variables, i);
                }
            }
        }
 
        if (my_reverse) {
            for (SeqLength i = 0; i < my_length; ++i) {
                char b = template_seq[my_length - i - 1];
                if (b != '-') {
                    add_base_to_hash(my_reverse_ref, complement_base(b));
                    add_mask_to_hash(my_reverse_mask);
                    my_reverse_mask_ambiguous.set(i);
                } else {
                    shift_hash(my_reverse_ref);
                    shift_hash(my_reverse_mask);
                    add_variable_base(my_reverse_variables, i);
                }
            }
        }
    }
 
public:
    struct State {
        SeqLength position = static_cast<SeqLength>(-1); // set to -1 so that the first call to next() overflows to 0.
 
        int forward_mismatches = 0;
 
        int reverse_mismatches = 0;
 
        bool finished = false;
 
        std::bitset<N> state;
        const char * seq;
        SeqLength len;
 
        std::bitset<N/4> ambiguous; // we only need a yes/no for the ambiguous state, so we can use a smaller bitset.
        SeqLength last_ambiguous; // contains the position of the most recent ambiguous base; should only be read if any_ambiguous = true.
        bool any_ambiguous = false; // indicates whether ambiguous.count() > 0.
    };
 
    State initialize(const char* read_seq, SeqLength read_length) const {
        State out;
        out.seq = read_seq;
        out.len = read_length;
 
        if (my_length <= read_length) {
            for (SeqLength i = 0; i < my_length - 1; ++i) {
                char base = read_seq[i];
 
                if (is_standard_base(base)) {
                    add_base_to_hash(out.state, base);
 
                    if (out.any_ambiguous) {
                        out.ambiguous <<= 1;
                    }
                } else {
                    add_other_to_hash(out.state);
 
                    if (out.any_ambiguous) {
                        out.ambiguous <<= 1;
                    } else {
                        out.any_ambiguous = true;
                    }
                    out.ambiguous.set(0);
                    out.last_ambiguous = i;
                }
            }
        } else {
            out.finished = true;
        }
 
        return out;
    }
 
    void next(State& state) const {
        SeqLength right = state.position + my_length;
        char base = state.seq[right];
 
        if (is_standard_base(base)) {
            add_base_to_hash(state.state, base); // no need to trim off the end, the mask will handle that.
            if (state.any_ambiguous) {
                state.ambiguous <<= 1;
 
                // If the last ambiguous position is equal to 'position', the
                // ensuing increment to the latter will shift it out of the
                // hash... at which point, we've got no ambiguity left.
                if (state.last_ambiguous == state.position) {
                    state.any_ambiguous = false;
                }
            }
 
        } else {
            add_other_to_hash(state.state);
 
            if (state.any_ambiguous) {
                state.ambiguous <<= 1;
            } else {
                state.any_ambiguous = true;
            }
            state.ambiguous.set(0);
            state.last_ambiguous = right;
        }
 
        ++state.position;
        full_match(state);
        if (right + 1 == state.len) {
            state.finished = true;
        }
 
        return;
    }
 
private:
    std::bitset<N> my_forward_ref, my_forward_mask;
    std::bitset<N> my_reverse_ref, my_reverse_mask;
    SeqLength my_length;
    bool my_forward, my_reverse;
 
    std::bitset<N/4> my_forward_mask_ambiguous; // we only need a yes/no for whether a position is an ambiguous base, so we can use a smaller bitset.
    std::bitset<N/4> my_reverse_mask_ambiguous;
 
    static void add_mask_to_hash(std::bitset<N>& current) {
        shift_hash(current);
        current.set(0);
        current.set(1);
        current.set(2);
        current.set(3);
        return;
    }
 
    static int strand_match(const State& match, const std::bitset<N>& ref, const std::bitset<N>& mask, const std::bitset<N/4>& mask_ambiguous) {
        // pop count here is equal to the number of non-ambiguous mismatches *
        // 2 + number of ambiguous mismatches * 3. This is because
        // non-ambiguous bases are encoded by 1 set bit per 4 bases (so 2 are
        // left after a XOR'd mismatch), while ambiguous mismatches are encoded
        // by all set bits per 4 bases (which means that 3 are left after XOR).
        int pcount = ((match.state & mask) ^ ref).count(); 
 
        if (match.any_ambiguous) {
            int acount = (match.ambiguous & mask_ambiguous).count();
            return (pcount - acount) / 2; // i.e., acount + (pcount - acount * 3) / 2;
        } else {
            return pcount / 2;
        }
    }
 
    void full_match(State& match) const {
        if (my_forward) {
            match.forward_mismatches = strand_match(match, my_forward_ref, my_forward_mask, my_forward_mask_ambiguous);
        }
        if (my_reverse) {
            match.reverse_mismatches = strand_match(match, my_reverse_ref, my_reverse_mask, my_reverse_mask_ambiguous);
        }
    }
 
private:
    std::vector<std::pair<SeqLength, SeqLength> > my_forward_variables, my_reverse_variables;
 
public:
    const std::vector<std::pair<SeqLength, SeqLength> >& forward_variable_regions() const {
        return my_forward_variables;
    }
 
    const std::vector<std::pair<SeqLength, SeqLength> >& reverse_variable_regions() const {
        return my_reverse_variables;
    } 
 
    const std::vector<std::pair<SeqLength, SeqLength> >& variable_regions(bool reverse) const {
        if (reverse) {
            return reverse_variable_regions();
        } else {
            return forward_variable_regions();
        }
    } 
 
private:
    static void add_variable_base(std::vector<std::pair<SeqLength, SeqLength> >& variables, SeqLength i) {
        if (!variables.empty()) {
            auto& last = variables.back().second;
            if (last == i) {
                ++last;
                return;
            }
        }
        variables.emplace_back(i, i + 1);
        return;
    }
};
 
}
 
#endif