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(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