src/LocARNA/multiple_alignment.hh

#ifndef LOCARNA_MULTIPLE_ALIGNMENT_HH
#define LOCARNA_MULTIPLE_ALIGNMENT_HH

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iosfwd>
#include <string>
#include <vector>
#include <map>

#include "aux.hh"
#include "string1.hh"
#include "scoring_fwd.hh"
#include "sequence_annotation.hh"

#include <assert.h>

#include <iostream>

namespace LocARNA {

    class Alignment;
    class AlignmentEdges;
    template <class T>
    class Alphabet;
    class BasePairs;
    class Scoring;
    class Sequence;

    class MultipleAlignment {
    public:
        typedef size_t size_type; 

        struct FormatType {
            enum type {
                STOCKHOLM, 
                PP,        
                CLUSTAL,   
                FASTA      
            };

            static size_t
            size() {
                return 4;
            }
        };

        struct AnnoType {
            enum type {
                consensus_structure,
                structure,
                fixed_structure,
                anchors
            };

            static size_t
            size() {
                return 4;
            }
        };

    private:
        typedef std::vector<std::vector<std::string> > annotation_tags_t;

        static annotation_tags_t annotation_tags;

        static void
        init_annotation_tags();

    public:
        static size_t
        num_of_annotypes() {
            return annotation_tags.size();
        }

        class SeqEntry {
        public:
            typedef MultipleAlignment::size_type size_type; 

            typedef std::pair<pos_type, pos_type>
                pos_pair_t; 

        private:
            std::string name_;        
            std::string description_; 
            string1 seq_;             //<! alignment string of the sequence

        public:
            SeqEntry(const std::string &name, const std::string &seq)
                : name_(name), description_(""), seq_((string1)seq) {}

            SeqEntry(const std::string &name, const string1 &seq)
                : name_(name), description_(""), seq_(seq) {}

            SeqEntry(const std::string &name,
                     const std::string &description,
                     const std::string &seq)
                : name_(name), description_(description), seq_((string1)seq) {}

            SeqEntry(const std::string &name,
                     const std::string &description,
                     const string1 &seq)
                : name_(name), description_(description), seq_(seq) {}

            // access

            const std::string &
            name() const {
                return name_;
            }

            const std::string &
            description() const {
                return description_;
            }

            const string1 &
            seq() const {
                return seq_;
            }

            size_type
            length_wogaps() const;

            //****************************************
            // projections

            pos_type
            pos_to_col(pos_type pos) const;

            pos_pair_t
            col_to_pos(pos_type col) const;

            void
            reverse() {
                seq_.reverse();
            }

            void
            push_back(char c) {
                seq_.push_back(c);
            }

            void
            set_seq(const string1 &seq) {
                seq_ = seq;
            }
        };

        class AliColumn {
            const MultipleAlignment &ma_;
            size_type col_index_;

        public:
            AliColumn(const MultipleAlignment &ma, size_type col_index)
                : ma_(ma), col_index_(col_index) {
                assert(1 <= col_index);
                assert(col_index <= ma.length());
            }

            const char &operator[](size_type row_index) const {
                return ma_.seqentry(row_index).seq()[col_index_];
            }

            size_type
            size() const {
                return ma_.num_of_rows();
            }

            bool
            operator==(const AliColumn &ac) const {
                bool ret = this->size() == ac.size();
                for (size_type i = 0; ret && i < size(); i++) {
                    ret = (this->ma_.seqentry(i).seq()[this->col_index_] ==
                           ac.ma_.seqentry(i).seq()[ac.col_index_]);
                }
                return ret;
            }

            bool
            operator!=(const AliColumn &ac) const {
                return !(*this == ac);
            }
        };

    private:
        typedef std::map<std::string, size_type> str2idx_map_t;

        typedef std::map<size_t, SequenceAnnotation> annotation_map_t;

        //************************************************************
        // attributes of MultipleAlignment

        std::vector<SeqEntry> alig_;

        annotation_map_t annotations_;

        str2idx_map_t name2idx_;

        // end attributes
        //************************************************************

        void
        create_name2idx_map();

        void
        read_clustallike(std::istream &in, FormatType::type format);

        void
        read_stockholm(std::istream &in);

        void
        read_clustalw(std::istream &in);

        void
        read_fasta(std::istream &in);

    public:
        typedef std::vector<SeqEntry>::const_iterator const_iterator;

        MultipleAlignment();

        MultipleAlignment(const std::string &file,
                          FormatType::type format = FormatType::CLUSTAL);

        MultipleAlignment(std::istream &in,
                          FormatType::type format = FormatType::CLUSTAL);

        MultipleAlignment(const std::string &name, const std::string &sequence);

        MultipleAlignment(const std::string &nameA,
                          const std::string &nameB,
                          const std::string &alistringA,
                          const std::string &alistringB);

        MultipleAlignment(const Alignment &alignment,
                          bool only_local = false,
                          bool special_gap_symbols = false);

        MultipleAlignment(const AlignmentEdges &edges,
                          const Sequence &seqA,
                          const Sequence &seqB);

    protected:
        void
        init(const AlignmentEdges &edges,
             const Sequence &seqA,
             const Sequence &seqB,
             bool special_gap_symbols);

    public:
        virtual ~MultipleAlignment();

        const Sequence &
        as_sequence() const;

        void
        normalize_rna_symbols();

        size_type
        num_of_rows() const {
            return alig_.size();
        }

        bool
        empty() const {
            return alig_.empty();
        }

        const SequenceAnnotation &
        annotation(const AnnoType::type &annotype) const;

        void
        set_annotation(const AnnoType::type &annotype,
                       const SequenceAnnotation &annotation) {
            assert(0 <= annotype && annotype < num_of_annotypes());
            annotations_[(size_t)annotype] = annotation;
        }

        bool
        has_annotation(const AnnoType::type &annotype) const {
            assert(0 <= annotype && annotype < num_of_annotypes());
            return annotations_.find(annotype) != annotations_.end();
        }

        bool
        is_proper() const;

        pos_type
        length() const {
            return alig_.empty() ? 0 : alig_[0].seq().length();
        }

        const_iterator
        begin() const {
            return alig_.begin();
        }

        const_iterator
        end() const {
            return alig_.end();
        }

        bool
        contains(const std::string &name) const;

        /* index access saves time over access by sequence name */

        size_type
        index(const std::string &name) const {
            str2idx_map_t::const_iterator it = name2idx_.find(name);
            assert(it != name2idx_.end());
            return it->second;
        }

        const SeqEntry &
        seqentry(size_type index) const {
            return alig_[index];
        }

        const SeqEntry &
        seqentry(const std::string &name) const {
            return alig_[index(name)];
        }

        size_type
        deviation(const MultipleAlignment &ma) const;

        double
        sps(const MultipleAlignment &ma, bool compalign = true) const;

        double
        cmfinder_realignment_score(const MultipleAlignment &ma) const;

        double
        avg_deviation_score(const MultipleAlignment &ma) const;

        std::string
        consensus_sequence() const;

        AliColumn
        column(size_type col_index) const {
            return AliColumn(*this, col_index);
        }

        void
        append(const SeqEntry &seqentry);

        void
        prepend(const SeqEntry &seqentry);

        void
        operator+=(const AliColumn &c);

        void
        operator+=(char c);

        void
        reverse();

        // ------------------------------------------------------------
        // output

        std::ostream &
        write(std::ostream &out,
              FormatType::type format =
                  MultipleAlignment::FormatType::CLUSTAL) const;

        std::ostream &
        write(std::ostream &out,
              size_t width,
              FormatType::type format =
                  MultipleAlignment::FormatType::CLUSTAL) const;

        std::ostream &
        write_name_sequence_line(std::ostream &out,
                                 const std::string &name,
                                 const std::string &sequence,
                                 size_t namewidth) const;

        std::ostream &
        write(std::ostream &out,
              size_type start,
              size_type end,
              FormatType::type format =
                  MultipleAlignment::FormatType::CLUSTAL) const;

        bool
        checkAlphabet(const Alphabet<char> &alphabet) const;

    private:
        static size_type
        deviation2(const string1 &a1,
                   const string1 &a2,
                   const string1 &ref1,
                   const string1 &ref2);

        static double
        pairwise_match_score(const SeqEntry &a1,
                             const SeqEntry &a2,
                             const SeqEntry &ref1,
                             const SeqEntry &ref2,
                             bool score_common_gaps);

        static std::vector<int>
        match_vector(const string1 &s, const string1 &t);

        static std::vector<int>
        match_vector2(const string1 &s, const string1 &t);

        static size_t
        count_matches(const SeqEntry &a1, const SeqEntry &a2);

        static size_t
        count_exclusive_matches(const SeqEntry &a1,
                                const SeqEntry &a2,
                                const SeqEntry &ref1,
                                const SeqEntry &ref2);

        static double
        pairwise_deviation_score(const SeqEntry &a1,
                                 const SeqEntry &a2,
                                 const SeqEntry &ref1,
                                 const SeqEntry &ref2);

    public:
        void
        write_debug(std::ostream &out = std::cout) const;
    };

    std::ostream &
    operator<<(std::ostream &out, const MultipleAlignment &ma);

} // end namespace

#endif // LOCARNA_MULTIPLE_ALIGNMENT_HH