SimpleSingleMatch.hpp

Go to the documentation of this file.

#ifndef KAORI_SIMPLE_SINGLE_MATCH_HPP
#define KAORI_SIMPLE_SINGLE_MATCH_HPP
 
#include "ScanTemplate.hpp"
#include "BarcodePool.hpp"
#include "BarcodeSearch.hpp"
#include "utils.hpp"
 
#include <string>
#include <unordered_map>
#include <vector>
 
namespace kaori {
 
template<SeqLength max_size_>
class SimpleSingleMatch {
public:
    struct Options {
        int max_mismatches = 0;
 
        DuplicateAction duplicates = DuplicateAction::ERROR;
 
        SearchStrand strand = SearchStrand::FORWARD;
    };
 
public:
    SimpleSingleMatch(
        const char* template_seq, 
        SeqLength template_length, 
        const BarcodePool& barcode_pool, 
        const Options& options
    ) : 
        my_forward(search_forward(options.strand)), 
        my_reverse(search_reverse(options.strand)),
        my_max_mm(options.max_mismatches),
        my_constant(template_seq, template_length, options.strand)
    {
        // Exact strandedness doesn't matter here, just need the number and length.
        const auto& regions = my_constant.forward_variable_regions();
        if (regions.size() != 1) {
            throw std::runtime_error("expected one variable region in the constant template");
        }
 
        SeqLength var_length = regions[0].second - regions[0].first;
        if (var_length != barcode_pool.length()) {
            throw std::runtime_error("length of barcode_pool sequences (" + std::to_string(barcode_pool.length()) + 
                ") should be the same as the barcode_pool region (" + std::to_string(var_length) + ")");
        }
 
        SimpleBarcodeSearch::Options bs_opt;
        bs_opt.duplicates = options.duplicates;
        bs_opt.max_mismatches = my_max_mm;
 
        if (my_forward) {
            bs_opt.reverse = false;
            my_forward_lib = SimpleBarcodeSearch(barcode_pool, bs_opt);
        }
        if (my_reverse) {
            bs_opt.reverse = true;
            my_reverse_lib = SimpleBarcodeSearch(barcode_pool, bs_opt);
        }
    }
 
private:
    bool my_forward, my_reverse;
    int my_max_mm;
    ScanTemplate<max_size_> my_constant;
    SimpleBarcodeSearch my_forward_lib, my_reverse_lib;
 
public:
    struct State {
        BarcodeIndex index = STATUS_UNMATCHED;
 
        SeqLength position = 0;
 
        int mismatches = 0;
 
        int variable_mismatches = 0;
 
        bool reverse = false;
 
        std::string buffer;
        typename SimpleBarcodeSearch::State forward_details, reverse_details;
    };
 
    State initialize() const {
        return State();
    }
 
    void reduce(State& state) {
        if (my_forward) {
            my_forward_lib.reduce(state.forward_details);
        }
        if (my_reverse) {
            my_reverse_lib.reduce(state.reverse_details);
        }
    }
 
private:
    void forward_match(const char* seq, const typename ScanTemplate<max_size_>::State& details, State& state) const {
        auto start = seq + details.position;
        const auto& range = my_constant.forward_variable_regions()[0];
        state.buffer.clear();
        state.buffer.insert(state.buffer.end(), start + range.first, start + range.second);
        my_forward_lib.search(state.buffer, state.forward_details, my_max_mm - details.forward_mismatches);
    }
 
    void reverse_match(const char* seq, const typename ScanTemplate<max_size_>::State& details, State& state) const {
        auto start = seq + details.position;
        const auto& range = my_constant.reverse_variable_regions()[0];
        state.buffer.clear();
        state.buffer.insert(state.buffer.end(), start + range.first, start + range.second);
        my_reverse_lib.search(state.buffer, state.reverse_details, my_max_mm - details.reverse_mismatches);
    }
 
public:
    bool search_first(const char* read_seq, SeqLength read_length, State& state) const {
        auto deets = my_constant.initialize(read_seq, read_length);
        bool found = false;
        state.index = STATUS_UNMATCHED;
        state.mismatches = 0;
        state.variable_mismatches = 0;
 
        auto update = [&](bool rev, int const_mismatches, const typename SimpleBarcodeSearch::State& x) -> bool {
            if (!is_barcode_index_ok(x.index)) {
                return false;
            }
 
            int total = const_mismatches + x.mismatches;
            if (total > my_max_mm) {
                return false;
            }
 
            found = true;
            state.position = deets.position;
            state.mismatches = total;
            state.reverse = rev;
            state.index = x.index;
            state.variable_mismatches = x.mismatches;
            return true;
        };
 
        while (!deets.finished) {
            my_constant.next(deets);
 
            if (my_forward && deets.forward_mismatches <= my_max_mm) {
                forward_match(read_seq, deets, state);
                if (update(false, deets.forward_mismatches, state.forward_details)) {
                    break;
                }
            }
 
            if (my_reverse && deets.reverse_mismatches <= my_max_mm) {
                reverse_match(read_seq, deets, state);
                if (update(true, deets.reverse_mismatches, state.reverse_details)) {
                    break;
                }
            }
        }
 
        return found;
    }
 
    bool search_best(const char* read_seq, SeqLength read_length, State& state) const {
        auto deets = my_constant.initialize(read_seq, read_length);
        state.index = STATUS_UNMATCHED;
        bool found = false;
        int best = my_max_mm + 1;
 
        auto update = [&](bool rev, int const_mismatches, const typename SimpleBarcodeSearch::State& x) -> void {
            if (!is_barcode_index_ok(x.index)) {
                return;
            }
 
            auto total = x.mismatches + const_mismatches;
            if (total == best) { 
                if (state.index != x.index) { // ambiguous, setting back to a mismatch.
                    found = false;
                    state.index = STATUS_AMBIGUOUS;
                }
            } else if (total < best) {
                found = true;
                best = total; 
                // A further optimization at this point would be to narrow
                // max_mm to the current 'best'. But this probably
                // isn't worth it.
 
                state.index = x.index;
                state.mismatches = total;
                state.variable_mismatches = x.mismatches;
                state.position = deets.position;
                state.reverse = rev;
            }
        };
 
        while (!deets.finished) {
            my_constant.next(deets);
 
            if (my_forward && deets.forward_mismatches <= my_max_mm) {
                forward_match(read_seq, deets, state);
                update(false, deets.forward_mismatches, state.forward_details);
            }
 
            if (my_reverse && deets.reverse_mismatches <= my_max_mm) {
                reverse_match(read_seq, deets, state);
                update(true, deets.reverse_mismatches, state.reverse_details);
            }
        }
 
        return found;
    }
};
 
}
 
#endif