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, size_t N>
    class Alphabet;
    class BasePairs;
    class Scoring;
    class Sequence;
 
    class MultipleAlignment {
    public:
        enum class FormatType {
            STOCKHOLM, 
            PP,        
            CLUSTAL,   
            FASTA      
        };
        static const std::vector<FormatType> FormatTypes;
 
        enum class AnnoType {
            consensus_structure,
            structure,
            fixed_structure,
            anchors
        };
        static const std::vector<AnnoType> AnnoTypes;
 
        static size_t
        num_of_annotypes() {
            return annotation_tags.size();
        }
 
        class SeqEntry {
        public:
 
            typedef std::pair<pos_type, pos_type>
                pos_pair_t; 
 
            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;
            }
 
        private:
            std::string name_;        
            std::string description_; 
            string1 seq_;             //<! alignment string of the sequence
 
        }; //end class SeqEntry
 
        class AliColumn {
        public:
            using value_type = char;
 
            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);
            }
 
            // make (forward) iterable
 
            class const_iterator {
            public:
 
                const_iterator
                operator ++() {
                    ++row_index_;
                    return *this;
                }
 
                const char &
                operator *() {
                    return col_[row_index_];
                }
 
                bool
                operator != (const const_iterator &it) {
                    return row_index_ != it.row_index_ || col_!=it.col_;
                }
 
            private:
                friend class AliColumn;
 
                const_iterator(const AliColumn &col, size_type row_index)
                    : col_(col), row_index_(row_index) {
                }
 
                const AliColumn &col_;
                size_type row_index_;
            };
 
            using iterator = const_iterator;
 
            auto
            begin() const {
                return const_iterator(*this,0);
            }
 
            auto
            end() const {
                return const_iterator(*this,this->size());
            }
 
        private:
            const MultipleAlignment &ma_;
            size_type col_index_;
        }; // end class AliColumn
 
    public:
        typedef std::vector<SeqEntry>::const_iterator const_iterator;
        typedef std::vector<SeqEntry>::iterator iterator;
 
        MultipleAlignment();
 
        MultipleAlignment(const MultipleAlignment &ma) = default;
 
        MultipleAlignment(MultipleAlignment &&ma) = default;
 
        MultipleAlignment &
        operator =(const MultipleAlignment &ma) = default;
        MultipleAlignment &
        operator =(MultipleAlignment &&ma) = default;
 
        virtual ~MultipleAlignment();
 
 
        MultipleAlignment(const std::string &file,
                          FormatType format = FormatType::CLUSTAL);
 
        MultipleAlignment(std::istream &in,
                          FormatType 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);
 
        void
        normalize_rna_symbols();
 
        size_type
        num_of_rows() const {
            return alig_.size();
        }
 
        bool
        empty() const {
            return alig_.empty();
        }
 
        const SequenceAnnotation &
        annotation(const AnnoType &annotype) const;
 
        void
        set_annotation(const AnnoType &annotype,
                       const SequenceAnnotation &annotation) {
            annotations_[annotype] = annotation;
        }
 
        bool
        has_annotation(const AnnoType &annotype) const {
            return annotations_.find(annotype) != annotations_.end();
        }
 
        bool
        is_proper() const;
 
        pos_type
        length() const {
            return alig_.empty() ? 0 : alig_[0].seq().length();
        }
 
        using value_type = SeqEntry;
 
 
        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 format =
                  MultipleAlignment::FormatType::CLUSTAL) const;
 
        std::ostream &
        write(std::ostream &out,
              size_t width,
              FormatType 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 format =
                  MultipleAlignment::FormatType::CLUSTAL) const;
 
        template <size_t N>
        bool
        checkAlphabet(const Alphabet<char, N> &alphabet) const;
 
        void
        write_debug(std::ostream &out = std::cout) const;
 
    protected:
        void
        init(const AlignmentEdges &edges,
             const Sequence &seqA,
             const Sequence &seqB,
             bool special_gap_symbols);
 
    private:
        class annotation_tags_t : public std::map< FormatType, std::map< AnnoType, std::string> > {
            //@brief constructor: initialize annotation tags table
        public:
            annotation_tags_t();
        };
 
        typedef std::map<std::string, size_type> str2idx_map_t;
 
        typedef std::map<AnnoType, SequenceAnnotation> annotation_map_t;
 
        //************************************************************
        // attributes of MultipleAlignment
 
        static annotation_tags_t annotation_tags;
 
        std::vector<SeqEntry> alig_;
 
        annotation_map_t annotations_;
 
        str2idx_map_t name2idx_;
 
        // end attributes
        //************************************************************
 
 
        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);
 
        void
        create_name2idx_map();
 
        void
        read_clustallike(std::istream &in, FormatType format);
 
        void
        read_stockholm(std::istream &in);
 
        void
        read_clustalw(std::istream &in);
 
        void
        read_fasta(std::istream &in);
    };
 
    std::ostream &
    operator<<(std::ostream &out, const MultipleAlignment &ma);
 
 
    template <size_t N>
    bool
    MultipleAlignment::checkAlphabet(const Alphabet<char, N> &a) const {
        return std::all_of(alig_.begin(), alig_.end(),
                           [&a](const auto &row) {
                               return std::all_of(row.seq().begin(), row.seq().end(),
                                                  [&a](const auto &c) {
                                                      return a.in(c);
                                                  });
                           });
    }
 
 
} // end namespace
 
#endif // LOCARNA_MULTIPLE_ALIGNMENT_HH