alignment.hh

/*
 * alignment.hh
 *
 * License: Artistic License, see file LICENSE.TXT or 
 *          https://opensource.org/licenses/artistic-license-1.0
 */
 
#ifndef _ALIGNMENT
#define _ALIGNMENT
 
#include "types.hh"
 
#include <vector>
#include <list>
#include <iostream>
#include <exception>
#include <climits>
 
class MsaSignature; // forward declaration
 
class fragment {
public:
    fragment(int chrPos, int aliPos, int len) : chrPos(chrPos), aliPos(aliPos), len(len) {}
    fragment() {}
    int chrEnd() const { return chrPos + len - 1; }
    int chrPos; // chromosomal start position of fragment, 0-based
    int aliPos; // start position of fragment in alignment, 0-based
    int len;    // fragment length
 
    /*
    *   added by Giovanna Migliorelli 14.05.2020 
    *   code responsible for serialization
    *   aknowledgement : https://www.boost.org/doc/libs/1_70_0/libs/serialization/doc/tutorial.html
    */
 
    #ifdef TESTING
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive & ar, const unsigned int version)
    {
        ar & chrPos; 
        ar & aliPos;
        ar & len; 
    }
    #endif    
};
 
class AlignmentRow {
public:
    AlignmentRow() : cumFragLen(0) {}
    AlignmentRow(string seqID, int chrPos, Strand strand, string rowbuffer);
    ~AlignmentRow(){}
 
    int chrStart() const;
    int chrEnd() const;
    int aliEnd() const;
    int getSeqLen() const { return chrEnd() - chrStart() + 1; }
    int getCumFragLen() const { return cumFragLen; }
    void setCumFragLen(int len) { cumFragLen = len;} // use with care to ensure consistency
    int gapLenAfterFrag(size_t i) const {
    if (i+1 >= frags.size())
        return 0;
    return frags[i+1].chrPos - frags[i].chrPos - frags[i].len;
    }
    void addFragment(int chrPos, int aliPos, int len);
    void addFragment(fragment &f) { addFragment(f.chrPos, f.aliPos, f.len); }
    string getSignature() const {return seqID + ((strand == minusstrand)? "-" : "+");}
    void pack();
    friend ostream& operator<< (ostream& strm, const AlignmentRow &row);
    friend void appendRow(AlignmentRow **r1, const AlignmentRow *r2, int aliLen1, string sigstr);
 
    int getAliPos(int chrPos, vector<fragment>::const_iterator from);
    int getAliPos(int chrPos, vector<fragment>::const_iterator *from); // variant from Patrick Balmerth
 
    int getAliPos(int chrPos) { return getAliPos(chrPos, frags.begin()); }
 
    // convert from alignment to chromosomal position (inverse function of getAliPos())
    int getChrPos(int aliPos, vector<fragment>::const_iterator from);
    int getChrPos(int aliPos) { return getChrPos(aliPos, frags.begin()); }
    
    // data members
    string seqID; // e.g. chr21
    
    #ifdef TESTING
    int seqIDarchive;   // GM temporarily added to come around some problem with string serialization
    int chrLen;         // GM added to properly set rsa chr lengths during deserialization
    #endif
 
    Strand strand;
    vector<fragment> frags; // fragments are sorted by alignment positions AND by chromosomal positions (assumption for now)
private:
    int cumFragLen; // total length of all fragments, becomes a nonredundant attribute after merging of alignments
 
    /*
    *   added by Giovanna Migliorelli 14.05.2020 
    *   code responsible for serialization
    *   Acknowledegment : https://www.boost.org/doc/libs/1_70_0/libs/serialization/doc/tutorial.html
    */
 
    #ifdef TESTING
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive & ar, const unsigned int version)
    {
        ar & cumFragLen;
        ar & seqIDarchive; 
        ar & strand;
        ar & frags; 
        ar & chrLen;    
    }
    #endif
};
 
 
 
class Alignment {
public:
    Alignment(size_t k) : aliLen(0), rows(k, NULL) {} // initialize with NULL, which stand for missing AlignmentRows
    ~Alignment(){
    for (unsigned i=0; i<rows.size(); i++) 
        delete rows[i]; 
    }
    friend bool mergeable (Alignment *a1, Alignment *a2, int maxGapLen, float mergeableFrac, bool strong);
    friend ostream& operator<< (ostream& strm, const Alignment &a);
    void printTextGraph(ostream& strm);
    void merge(Alignment *other, const MsaSignature *sig = NULL); // append 'other' Alignment to this
    friend Alignment* mergeAliList(list<Alignment*> alis,  const MsaSignature *sig);
    friend void capAliSize(list<Alignment*> &alis, int maxRange);
    friend void reduceOvlpRanges(list<Alignment*> &alis, size_t maxCov, float covPen);
    int maxRange(); // chromosomal range, maximized over rows
    int numRows() const { return rows.size(); }
    int numFilledRows() const; // number of nonempty rows
    int getCumFragLen() const; // total length of all fragments
    int getMaxSeqIdLen() const;
    friend int medianChrStartEndDiff(Alignment *a1, Alignment *a2);
    string getSignature() const;
    int numFitting(const MsaSignature *sig) const;
    void shiftAliPositions(int offset);
    void pack(); // merge pairs of fragments without gap into one fragment making the alignment representation more compact
 
    #ifdef TESTING
    void convertAlignment(int r, int start, int end);
    #endif
 
public: // should rather be private
    int aliLen; // all aligned sequences are this long when gaps are included
    vector<AlignmentRow*> rows;
 
private:
    /*
    *   added by Giovanna Migliorelli 14.05.2020 
    *   code responsible for serialization and default constructor required in serialization
    *   aknowledgement : https://www.boost.org/doc/libs/1_70_0/libs/serialization/doc/tutorial.html
    */
 
    #ifdef TESTING
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive & ar, const unsigned int version)
    {
        ar & aliLen;
        ar & rows;
    }
    
    Alignment(){};
    #endif
};
 
int medianChrStartEndDiff(Alignment *a1, Alignment *a2);
 
// sorting operator, with respect to a given species index
struct SortCriterion {
    SortCriterion(size_t speciesIdx) : s(speciesIdx) {};
    bool operator() (Alignment* const a1, Alignment* const a2){
    // alignments, where row s is missing come last
    if (a2->rows[s] == NULL)
        return true;
    if (a1->rows[s] == NULL)
        return false;
    if (a1->rows[s]->seqID < a2->rows[s]->seqID)
        return true;
    if (a1->rows[s]->seqID > a2->rows[s]->seqID)
        return false;
    // same sequence, compare chromosomal start positions
    return (a1->rows[s]->chrStart() < a2->rows[s]->chrStart());
    }
    size_t s;
};
 
 
// sorting operator, with respect to a given species index, index version
// for sorting a vector or list of indices to another vector that holds the Alignments
struct IdxSortCriterion {
    IdxSortCriterion(vector<Alignment*> const &a_, size_t speciesIdx) : s(speciesIdx), a(a_) {};
    bool operator() (int i, int j){
    // alignments, where row s is missing come last
    if (a[j]->rows[s] == NULL && a[i]->rows[s]) // this conjunction makes sorting stable where the sequence is missing
        return true;
    if (a[i]->rows[s] == NULL)
        return false;
    if (a[i]->rows[s]->seqID < a[j]->rows[s]->seqID)
        return true;
    if (a[i]->rows[s]->seqID > a[j]->rows[s]->seqID)
        return false;
    // same sequence, compare chromosomal start positions
    return (a[i]->rows[s]->chrStart() < a[j]->rows[s]->chrStart());
    }
    size_t s;
    vector<Alignment*> const &a;
};
 
/*
 * b1 and b2 can be merged in that order because they are very similar and right next to each other.
 * In at least 'mergeableFrac' of the alignment rows the aligned sequenes are
 * strong=false: - refer to the same terget sequence, are on the same strand and satisfy 0 <= gaplen <= maxGapLen
 *                 if both are present
 * strong=true:  - refer to the same terget sequence, are on the same strand and satisfy 0 <= gaplen <= maxGapLen
 *                 if at least one is present
 */
bool mergeable (Alignment *b1, Alignment *b2, int maxGapLen, float mergeableFrac, bool strong);
 
inline bool isGap(char c){
    return (c == '-' || c == '.');
}
 
class MsaSignature {
public:
    string sigstr() const{
    string str;
    for (unsigned s = 0; s < sigrows.size(); ++s)
        if (sigrows[s] != "")
        str += itoa(s) + ":" + sigrows[s];
    return str;
    }
    vector<string> sigrows; // each row contains seqID and strand, e.g. chr21+
    int numAli;
    int sumAliLen;
    int sumCumFragLen;
    int depth;
    int color;
    bool operator< (const MsaSignature &other) const {
    return (sumCumFragLen > other.sumCumFragLen);
    // before: sorting by 1) number of species and 2) sumAliLen
    // return (depth > other.depth || (depth == other.depth && sumAliLen > other.sumAliLen));
    }
    bool fits(const Alignment &a, size_t s) const {
    return a.rows[s] && (sigrows[s] == a.rows[s]->seqID + strandChar(a.rows[s]->strand));
    }
    static int maxSigStrLen;
};
 
bool cmpSigPtr(const MsaSignature *s1, const MsaSignature *s2);
 
typedef pair<size_t, int> BoundaryFragment; // (s, chrPos), where s= species index
 
class CompareBoundaryFragment {
public:
    bool operator()(BoundaryFragment& bf1, BoundaryFragment& bf2) {
    return bf1.second > bf2.second; // => sort by increasing chromosomal position
    }
};
 
 
class MsaInsertion {
public:
    MsaInsertion(size_t s, size_t insertpos, string insert): s(s), insertpos(insertpos), insert(insert){ }
    MsaInsertion() { }
    bool operator<(const MsaInsertion& other) const {
        if (insertpos > other.insertpos)
            return true;
        if (insertpos < other.insertpos)
            return false;
        if (insert.length() > other.insert.length())
            return true;
        if (insert.length() < other.insert.length())
            return false;
        if (s < other.s)
            return true;
        return false;
    }
    size_t s; 
    size_t insertpos; 
    string insert; 
};
 
 
class StringAlignment {
public:
    StringAlignment(size_t numrows) : rows(numrows, ""), k(numrows), len(0) {} // initialize with empty strings
    ~StringAlignment(){len = 0;}
 
    void insert(std::list<MsaInsertion> &insList, int maxInsertLen = INT_MAX);
 
    size_t removeGapOnlyCols();
 
    bool isGapOnlyCol(size_t col){
        for (size_t s = 0; s < k; ++s)
            if (!rows[s].empty() && rows[s].at(col) != '-')
                return false;
        return true;
    }
 
    void computeLen(){
        size_t m = 0, rowlen;
        for (size_t s = 0; s < k; ++s){
            rowlen = rows[s].length();
            if (rowlen > 0){
                if (m == 0)
                    m = rowlen;
                else if (m != rowlen)
                    throw length_error("StringAlignment with rows of differing lengths");
            }
        }
        len = m;
    }
 
    friend ostream& operator<< (ostream& strm, const StringAlignment &msa);
 
// data members
    vector<string> rows;
    size_t k;
    size_t len;
};
 
 
#endif  // _ALIGNMENT