align.hpp

Go to the documentation of this file.

/*
 * Copyright 2008,2009,2010,2011,2012 Stephan Ossowski, Korbinian Schneeberger,
 * Felix Ott, Joerg Hagmann, Alf Scotland, Sebastian Bender
 * 
 * This file is part of SHORE.
 * 
 * SHORE is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * SHORE is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with SHORE.  If not, see <http://www.gnu.org/licenses/>.
 */


 
#ifndef SHORE_ALGO_ALIGN_HPP__
#define SHORE_ALGO_ALIGN_HPP__

#include <iterator>
#include <string>
#include <vector>

#include "shore/datatype/alignment_string_ops.hpp"
#include "shore/datatype/datatypes.hpp"
#include "shore/datatype/dna_iterator.hpp"
#include "shore/stream/streams.hpp"
#include "shore/processing/pipeline.hpp"
#include "shore/processing/pipe_box.hpp"
#include "shore/processing/buffer_chain.hpp"
#include "shore/processing/filter.hpp" // conversion
#include "shore/processing/feed.hpp"
#include "shore/processing/extractor.hpp"

#include "shore/base/util.hpp"
#include "shore/base/stringops.hpp"

namespace shore {

enum AlignmentMode
{
        ALN_NONE=-1,
        ALN_GLOBAL=0,
        ALN_DANGLING_REF=1,
        ALN_DANGLING_QRY=2,
        ALN_DANGLING_ANY=3,
        ALN_LOCAL=7
};

S_MAKE_ENUM_TRAITS(AlignmentMode)


class dp_nuc_scoring
{
 private:

        int m_scoringmatrix[256];

 public:

        dp_nuc_scoring();
        dp_nuc_scoring(const dp_nuc_scoring &other);
        dp_nuc_scoring & operator=(const dp_nuc_scoring &other);

        int get_pair_score(const shore::nuc::base ref,
                           const shore::nuc::base qry) const;

        int get_ref_score(const shore::nuc::base ref) const;

        int get_qry_score(const shore::nuc::base qry) const;


        void set_scores(const int match,const int mismatch,const int gap);
        void set_scores(const int match,const int mismatch,const int gap,const int ambmatch);
        void read_scores(const std::string &fn,const int defaultscore=-1);

        int get_best_score() const;

        void print(std::ostream &out) const;
};


class dp_alignment_matrix
{
 private:

        std::vector<int> m_matrix;

        size_t m_size;

        int m_ncol;

        int m_nrow;

        int m_globalmax;
        int m_lastrowmax;
        int m_lastcolmax;

        std::vector<size_t> m_globalmax_pos;
        std::vector<size_t> m_lastrowmax_pos;
        std::vector<size_t> m_lastcolmax_pos;

 public:

        class const_iterator
        {
         protected:

                friend class dp_alignment_matrix;

                const dp_alignment_matrix *m_const_that;

                int m_index;

                const_iterator(const dp_alignment_matrix *matrix);

         public:

                const_iterator();

                const_iterator(const const_iterator &other);

                const_iterator & operator=(const const_iterator &other);

                int get() const;

                int get_prev() const;

                int get_ref_prev() const;

                int get_qry_prev() const;

                void next();

                void ref_next();

                void qry_next();

                void ref_return();

                void qry_return();

                void prev();

                void ref_prev();

                void qry_prev();

                bool is_begin() const;

                bool is_end() const;

                size_t get_ref_pos() const;

                size_t get_qry_pos() const;

                size_t get_field_index() const;

                template<typename RefValue,typename QryValue,typename Scoring>
                int get_best_score(const RefValue & rv,const QryValue & qv,
                                   const Scoring & scoring,
                                   const bool local);
        };

        class iterator
        :public const_iterator
        {
         private:

                friend class dp_alignment_matrix;

                dp_alignment_matrix *m_that;


                iterator(dp_alignment_matrix *matrix);

         public:

                iterator();

                iterator(const iterator &other);

                iterator & operator=(const iterator &other);

                void set(const int value);
                void set_ref_last(const int value);
                void set_qry_last(const int value);
                void set_last(const int value);
        };

        dp_alignment_matrix();

        void reset_maxima();

        void resize(const size_t ref_size,const size_t qry_size);

        int get_score(const AlignmentMode right_ends) const;

        size_t num_global_max() const;

        size_t num_ref_max() const;

        size_t num_qry_max() const;

        int global_max() const;

        int ref_max() const;

        int qry_max() const;

        iterator begin();

        const_iterator last() const;

        const_iterator global_max_pos(size_t index) const;

        const_iterator ref_max_pos(size_t index) const;

        const_iterator qry_max_pos(size_t index) const;

        template<typename RefIter,typename QryIter>
        void print(std::ostream &os,
                   RefIter ir,QryIter iq) const;
};


template<typename RefType,typename QryType>
struct dp_traits
{
        typedef void scoring_type;
        typedef void alignment_builder_type;
};

template<>
struct dp_traits<shore::nuc::base,shore::nuc::base>
{
        typedef dp_nuc_scoring scoring_type;
        typedef shore::alignment_builder alignment_builder_type;
};


template<typename RefIter,typename QryIter,typename Scoring>
void dp_align(dp_alignment_matrix & matrix,
              RefIter rbeg,const size_t rsize,
              QryIter qbeg,const size_t qsize,
              const Scoring scoring,const AlignmentMode left_ends,
              const AlignmentMode exclude_left=ALN_NONE,
              const AlignmentMode exclude_right=ALN_NONE);


template<typename SeqIter>
struct dp_sequence
{
        SeqIter beg;
        SeqIter end;
        size_t ofs;
        size_t size;

        dp_sequence()
        :ofs(0),size(0)
        {}

        dp_sequence(const SeqIter & b,
                    const SeqIter & e,
                    const int of=0,const size_t sz=0)
        :beg(b),end(e),ofs(of),size(sz)
        {
                if(size==0)
                        size=e-b;
                else
                        size=std::min(size_t(e-b),size);
        }
};


template<typename RefIter,typename QryIter=RefIter,
         typename Scoring=typename dp_traits<typename std::iterator_traits<RefIter>::value_type,
                                             typename std::iterator_traits<QryIter>::value_type>::scoring_type>
struct dp_alignment
{
        AlignmentMode left_ends;

        AlignmentMode right_ends;

        const Scoring *scoring;

        dp_sequence<RefIter> ref;

        dp_sequence<QryIter> qry;

        int ref_id;

        int ref_hits;

        dp_alignment_matrix alignment_matrix;


        dp_alignment()
        :left_ends(ALN_DANGLING_QRY),
         right_ends(ALN_LOCAL),
         scoring(0),
         ref_id(0),
         ref_hits(0)
        {}
};


struct dp_aligner_config
{
        bool report_all_best;

        std::vector<AlignmentMode> left_ends;

        std::vector<AlignmentMode> right_ends;


        dp_aligner_config()
        :report_all_best(false)
        {}
};

template<typename RefIter,typename QryIter=RefIter,
         typename Scoring=typename dp_traits<typename std::iterator_traits<RefIter>::value_type,
                                             typename std::iterator_traits<QryIter>::value_type>::scoring_type>
class dp_aligner_pipe
:public shore::pipe_facade<dp_aligner_pipe<RefIter,QryIter,Scoring>,
                           dp_sequence<QryIter>,
                           dp_alignment<RefIter,QryIter,Scoring> >,
 public buffer_chain<dp_alignment<RefIter,QryIter,Scoring> >
{
 public:

        typedef dp_sequence<QryIter> append_type;
        typedef dp_alignment<RefIter,QryIter,Scoring> current_type;

        typedef buffer_chain<current_type> buffer_chain_type;

 private:

        friend class pipeline_core_access;

        dp_aligner_config m_conf;

        AlignmentMode m_end_modes[8];

        Scoring m_scoring;

        dp_sequence<RefIter> m_ref;

        dp_sequence<QryIter> m_qry;

        dp_alignment_matrix *m_alignment_matrix;

        int m_score;

        AlignmentMode m_left_ends;
        AlignmentMode m_right_ends;

        dp_alignment<RefIter,QryIter,Scoring> *m_auxiliary[8];

        std::vector<AlignmentMode> m_active_modes;


        void add_mode(const AlignmentMode left,const AlignmentMode right);


        bool has_weaker_constraint(const AlignmentMode left_a,
                                   const AlignmentMode left_b) const;

        void add_to_chain();

        void run();

        void next();

        void append(const dp_sequence<QryIter> & query);

        void flush();

 public:

        dp_aligner_pipe(const dp_aligner_config & c);

        dp_aligner_pipe(const dp_aligner_pipe & other);

        ~dp_aligner_pipe();

        void set_ref(RefIter rb,RefIter re);

        const dp_sequence<RefIter> get_ref() const;

        const Scoring & get_scoring() const;

        void set_scoring(const Scoring & matrix);
};


template<typename RefIter,typename QryIter=RefIter,
         typename Scoring=typename dp_traits<typename std::iterator_traits<RefIter>::value_type,
                                             typename std::iterator_traits<QryIter>::value_type>::scoring_type>
class dp_multialigner_pipe
:public shore::pipe_facade<dp_multialigner_pipe<RefIter,QryIter,Scoring>,
                           dp_sequence<QryIter>,
                           dp_alignment<RefIter,QryIter,Scoring> >,
 public buffer_chain<dp_alignment<RefIter,QryIter,Scoring> >
{
 public:

        typedef dp_sequence<QryIter> append_type;
        typedef dp_alignment<RefIter,QryIter,Scoring> current_type;

        typedef buffer_chain<current_type> buffer_chain_type;

 private:

        friend class pipeline_core_access;

        struct alignment_grabber
        :public shore::plugin<current_type>
        {
                dp_multialigner_pipe<RefIter,QryIter,Scoring> *that;
                int id;

                alignment_grabber(dp_multialigner_pipe<RefIter,QryIter,Scoring> *t,
                                  int i)
                :that(t),id(i)
                {}

                virtual int apply(std::deque<current_type *> & used,
                                  std::vector<current_type *> & unused,
                                  bool flush)
                {
                        if(used.empty())
                                return shore::plugin<current_type>::PLUGIN_STARVED;

                        const int score=
                                used.front()->alignment_matrix.get_score(used.front()->right_ends);

                        // Steal a buffer.
                        if(score>=that->m_score)
                        {
                                if(score>that->m_score)
                                {
                                        that->m_unused.insert(that->m_unused.end(),
                                                              that->m_used.begin(),
                                                              that->m_used.end());
                                        that->m_used.clear();
                                        that->m_score=score;
                                }
                                if(that->m_conf.report_all_best
                                   ||that->m_used.empty())
                                {
                                        used.front()->ref_id=id;
                                        that->m_used.push_back(used.front());
                                }
                                else
                                        that->m_unused.push_back(used.front());
                        }
                        else
                                that->m_unused.push_back(used.front());

                        used.pop_front();

                        // Return a buffer.
                        if(!that->m_unused.empty())
                        {
                                unused.push_back(that->m_unused.back());
                                that->m_unused.pop_back();
                        }

                        return shore::plugin<current_type>::PLUGIN_SUCCESS;
                }
        };

        dp_aligner_config m_conf;

        Scoring m_scoring;

        std::vector<dp_aligner_pipe<RefIter,QryIter,Scoring> *> m_aligners;

        std::deque<current_type *> m_used;
        std::vector<current_type *> m_unused;

        std::vector<alignment_grabber *> m_grabbers;

        int m_score;
        size_t m_ref_hits;

        shore::feed<dp_sequence<QryIter> > m_feed;


        void push_activated();

        void next();

        void append(const dp_sequence<QryIter> & query);

        void flush();

 public:

        dp_multialigner_pipe(const dp_aligner_config & c);

        dp_multialigner_pipe(const dp_multialigner_pipe & other);

        ~dp_multialigner_pipe();

        void add_ref(RefIter rb,RefIter re,const int id);

        const Scoring & get_scoring() const;

        void set_scoring(const Scoring & matrix);
};


struct dp_trace
{
        int alignment_id;

        AlignmentMode left_ends;

        AlignmentMode right_ends;

        std::string alignment_string;

        int score;

        int ref_id;

        int ref_hits;

        size_t ref_totalsize;

        size_t qry_totalsize;

        size_t ref_pos;

        size_t ref_end;

        size_t qry_pos;

        size_t qry_end;

        bool soft_clipped;

        bool unique_alignment;

        ExtendedFlag unique_mapping;


        dp_trace()
        :alignment_id(0),
         left_ends(ALN_DANGLING_QRY),
         right_ends(ALN_LOCAL),
         score(0),
         ref_id(0),
         ref_hits(0),
         ref_totalsize(0),
         qry_totalsize(0),
         ref_pos(0),ref_end(0),
         qry_pos(0),qry_end(0),
         soft_clipped(false),
         unique_alignment(false),
         unique_mapping(XFLAG_NA)
        {}
};

struct dp_backtracer_config
{
        bool trace_all_alignments;

        bool trace_all_mappings;

        bool check_unique_mapping;

        bool soft_clip;

        bool build_alignment;


        dp_backtracer_config()
        :trace_all_alignments(false),
         trace_all_mappings(false),
         check_unique_mapping(false),
         soft_clip(false),
         build_alignment(true)
        {}
};

template<typename RefIter,typename QryIter=RefIter,
         typename Scoring=typename dp_traits<typename std::iterator_traits<RefIter>::value_type,
                                             typename std::iterator_traits<QryIter>::value_type>::scoring_type,
         typename AlnBuilder=typename dp_traits<typename std::iterator_traits<RefIter>::value_type,
                                                typename std::iterator_traits<QryIter>::value_type>::alignment_builder_type>
class dp_backtracer_pipe
:public pipe_facade<dp_backtracer_pipe<RefIter,QryIter,Scoring,AlnBuilder>,
                    dp_alignment<RefIter,QryIter,Scoring>,
                    dp_trace>,
 public buffer_chain<dp_trace>
{
 private:

        friend class pipeline_core_access;

        enum TraceOp
        {
                TRACE_MATCH=1,
                TRACE_REFGAP=2,
                TRACE_QRYGAP=4
        };

        struct trace_state
        {
                dp_alignment_matrix::const_iterator iter;

                TraceOp trace_op;

                trace_state(const dp_alignment_matrix::const_iterator &it,
                            const TraceOp tr)
                :iter(it),trace_op(tr)
                {}
        };

        dp_backtracer_config m_conf;

        const dp_alignment<RefIter,QryIter,Scoring> *m_alignment;

        int m_alignment_id;

        size_t m_trace_index;

        size_t m_trace_size;

        size_t m_left_count;

        int m_trace_score;
        int m_trace_ref_end;
        int m_trace_qry_end;
        bool m_trace_aln_unique;
        ExtendedFlag m_trace_map_unique;

        AlnBuilder m_ws;

        std::vector<trace_state> m_trace_stack;
        std::vector<trace_state> m_trace_stack_save;

        bool m_assess_leftends_mode;

        std::vector<bool> m_visited;
        std::vector<bool> m_visited_save;

        bool m_trace_all_branches;

        void set_trace_size(const size_t size);

        ExtendedFlag mark_visited(const dp_alignment_matrix::const_iterator &it);

        void init_trace(const dp_alignment_matrix::const_iterator &it);

        // \brief  Start a backtrace from the next end position.
        // \return     Whether a new trace is ready to emit.
        bool start_trace();

        // \brief  Collect all alignments, put alternative paths on the stack.
        // \return     Whether a new trace is ready to emit.
        bool continue_trace();

        void assess_leftends();


        void next();

        void append(const dp_alignment<RefIter,QryIter,Scoring> & a);

        void flush();

 public:

        dp_backtracer_pipe(const dp_backtracer_config & c);
};


class dp_trace_filter
:public plugin<dp_trace>
{
 public:

        struct config
        {
                double threshold_rel;
                double threshold_abs;
                bool threshold_global;
                bool reject_ambiguous;

                config();
        };

 private:

        friend class pipeline_core_access;

        config m_conf;

        int m_best_score;

        virtual int apply(std::deque<dp_trace *> & used,
                          std::vector<dp_trace *> & unused,
                          const bool flushed);

 public:

        dp_trace_filter(const int best_score,const config & c);
};

struct dp_mapping
{
        const shore::read *read;
        const dp_trace *trace;

        dp_mapping()
        :read(0),
         trace(0)
        {}
};

class dp_mapper_pipe
:public shore::pipe_facade<dp_mapper_pipe,shore::read,dp_mapping>
{
 public:

        struct config
        {
                dp_aligner_config align_conf;

                dp_backtracer_config trace_conf;

                dp_trace_filter::config filter_conf;
        };

 private:

        friend class shore::pipeline_core_access;

        config m_conf;

        dp_nuc_scoring m_scoring;

        shore::feed<dp_sequence<shore::dna_iterator> > m_feed;
        dp_multialigner_pipe<shore::dna_iterator> m_aligner;
        dp_backtracer_pipe<shore::dna_iterator> m_tracer;
        dp_trace_filter m_filter;
        shore::extractor<dp_trace> m_extractor;

        int m_ref_id;

        dp_mapping m_current;

        void append(const shore::read &r);
        void next();

 public:

        dp_mapper_pipe(const int best_score,const config &c);
        dp_mapper_pipe(const dp_mapper_pipe & other);
        ~dp_mapper_pipe();

        void add_ref(shore::dna_iterator rb,shore::dna_iterator re);
        void add_ref(shore::dna_iterator rb_fwd,shore::dna_iterator re_fwd,
                     shore::dna_iterator rb_rev,shore::dna_iterator re_rev);

        void set_scoring(const dp_nuc_scoring & matrix);
};

typedef shore::pipe_box<dp_mapper_pipe> dp_mapper;


class dp_mapping_writer
{
 private:

        shore::ostreams m_os;

 public:

        typedef dp_mapping append_type;

        dp_mapping_writer(const std::string & outfile);

        void append(const dp_mapping & m);
        void flush();
};

typedef shore::sink<dp_mapping_writer> dp_mapping_sink;

} // namespace shore

// Include inline and template implementation.
#include "align.i.hpp"
#include "align.t.hpp"

#endif // SHORE_ALGO_ALIGN_HPP__