dynaln.h

Go to the documentation of this file.

#ifndef DYNALN_H
#define DYNALN_H

#include "bioseq.h"
#include "matrix.h"
#include <list>
#include <vector>

/* experimental, might not be needed 
struct Alignedcol {
   Alignedcol() { }
   Alignedcol(const char t, const char m, const char b) 
      : top(t), middle(m), buttom(b) { }
   char top;
   char middle;
   char buttom;
};

struct Alnstring {
   string top;
   string middle;
   string buttom;
};

struct Alnsummary {
   int score;
   int idencnt;
   int simcnt;  // excluding identical residues
   int numgaps1;
   int numgaps2;
   int gaplen1;
   int gaplen2;
};
*/

enum ALNTYPE {GLOBAL, LOCAL};
//enum POINTER {diag, top, left}; // for traceback

class Dynaln {
   public:
      Dynaln() : ST(), gapi(-8), gape(-2), seq1(0),
                 seq2(0),alnidx(), S(0), M(0),
                 IX(0), IY(0), Smax(0), idencnt(0), simcnt(0), numgaps1(0),
                 numgaps2(0), gaplen1(0), gaplen2(0),
                 seq1begin(-1), seq1end(-1), seq2begin(-1), seq2end(-1),
                 topaln(), middle(), buttomaln(), topruler(), buttomruler() { }
      Dynaln(const bioseq &s1, const bioseq &s2) 
         : ST(), gapi(-8), gape(-2), seq1(&s1), seq2(&s2), 
            alnidx(), S(0), M(0), IX(0), IY(0),
            Smax(0), idencnt(0), simcnt(0), numgaps1(0),
            numgaps2(0), gaplen1(0), gaplen2(0),
            seq1begin(-1), seq1end(-1), seq2begin(-1), seq2end(-1),
            topaln(), middle(), buttomaln(), topruler(),buttomruler()
            { if ( (s1.getSequenceType() & NUCLEIC_ACID & s2.getSequenceType()) > 0)
               ST.setMatrix("NUC.4.4");
                  }
      Dynaln(const bioseq &s1, const bioseq &s2, const string &ma) 
         : ST(ma), seq1(&s1), seq2(&s2), 
            alnidx(), S(0), M(0), IX(0), IY(0),
            Smax(0), idencnt(0), simcnt(0), numgaps1(0),
            numgaps2(0), gaplen1(0), gaplen2(0),
            seq1begin(-1), seq1end(-1), seq2begin(-1), seq2end(-1),
            topaln(), middle(), buttomaln(), topruler(), buttomruler()
            { 
               gapi=2*ST.getMinScore(); 
               gape=ST.getMinScore()/3; 
               if (gape > -1) gape = -1; 
               if (gape < -3) gape = -3; }
      int global();
      int runglobal() {
         int s=global(); buildResult(); return s; }
      int local();
      int runlocal() {
         int s=local(); buildResult(); return s; }
      //void globalTraceBack();
      ~Dynaln();

      void buildResult();
      void buildAlnInfo();


      void printAlign(ostream &ous, const int w=70) const;
      const string& getTopAln() const { return topaln; }
      const string& getButtomAln() const { return buttomaln; }
      void getNgMatchArray(vector<int> &ngmarr, bool append=false, bool countsim=true) const;
      void getMatchArray(vector<int> &marr, bool countsim=true) const;
      string toDelimitedString(const string &dl=",", int ibase=0) const;

      void showseq() const {
         cout << seq1->toString() << "\n"
            << seq2->toString() << endl; }

      int getScore() { return Smax; }
      int getAlnlen() const { return alnidx.size(); }
      float getIdentity() const { return (float)idencnt/getAlnlen(); }
      float getSimilarity() const { return (float)(idencnt+simcnt)/getAlnlen(); }
      pair<int,int> numgaps() const { 
         return make_pair(numgaps1, numgaps2); }
      static string headers(const string &dl);
      //const Matrix& getScoringMatrix() const { return ST; }
      const Matrix* getScoringMatrix() const { return &ST; }
      int topBeginIndex() const { return seq1begin; }
      int buttomBeginIndex() const { return seq2begin; }
      int topEndIndex() const { return seq1end; }
      int buttomEndIndex() const { return seq2end; }

      void setMatrix(const string &mat) { ST.setMatrix(mat); }
      void setMatrix(const Matrix &mat) { ST=mat; }

      void setSeq1(const bioseq &sq) { seq1 = &sq;  }
      void setSeq2(const bioseq &sq) { seq2 = &sq; }
      void setSeq(const bioseq &sq1, const bioseq &sq2) {
         seq1=&sq1; seq2=&sq2; }
      void setseq(const bioseq &sq1, const bioseq &sq2) {
         seq1=&sq1; seq2=&sq2; }
      //void setSeq(const string &sq1, const string &sq2) {
       //  *seq1=sq1; *seq2=sq2; }

      void setGapInsert(int ins) { gapi=ins; }
      void setGapExtend(int ext) { gape=ext; }
      void setGapParameter(const int go, const int ge) {
         gapi=go; gape=ge; }
      void setAllParameter(const Matrix& mat, const int go, const int ge)
      { ST=mat; gapi=go; gape=ge; }

      static const char gapchar='-';
      static const char idenchar='|';
      static const char simchar=':';
      static const int simcut=0;
      static const int PTRNULL=0;
      static const int PTRDIAG=1;
      static const int PTRTOP=2; // Ix
      static const int PTRLEFT=4; // Iy

      void debug_showmatrix(ostream &ous) const;

   protected:
      void countnumgaps();
      void clearResult();
      //void allocmemLS();
      
      /* Nr is the number of row (seq1 length + 1)
       * Nc is the number of columns (seq2 length + 1)
       */
      //void allocmem(const int Nr, const int Nc);
      void allocmem();
      //void traceback(int &i, int &j, int *&Sc, bool local=false);
      
      /*
       * the traceback function produce the alignment result
       * and put it into the alnidx list structure.
       * It also set the begin,end of each sequence that.
       * For global alignment, it will traceback the begining
       * gaps. The ending gap will also be removed from the
       * end of the sequence.
       *
       * i, and j are the matrix location index 0-based
       */
      //void traceback(int &i, int &j, bool local=false);
      //void traceback(int &i, int &j);
      void tracepointer(int &i, int &j);
      void findAlnBoundary();
      //void initialize();
      int gapi, gape;  // should be negative numbers
      //Matrix ST; // scoring table, use protein BLOSUM50 as default.
      //We should use pointers
      Matrix ST; // scoring table, use protein BLOSUM50 as default.
      // for DNA use NUC.4.4
      const bioseq *seq1;
      const bioseq *seq2;
      // the code were assigned when allocating memory
      const int* C1; // needs to be updated properly when sequence got updated
      const int* C2; // code for the sequence
      int* S;    // simulate 2-D array, 
                 //size (seq1.length()+1)*(seq2.length()+1)
      int Ssize; // for the next round, memorize the last state
      int *M, *IX, *IY; // array for computing the scores
      //int* Itop;  // for temporary usage
      int numcol;  // for the next round
      //int* P; // for trace back pointer
      //int Ileft;  // for temporary usage
      //stack<pair<int, int> > alnidx;
      
      list<pair<int, int> > alnidx;
      //int score;  // assigned in the trace back step.
      int Smax;  // assigned before the trace back step.
      // in local it is use during the build up phase
      // in global it is the last cell m,n
      int Smaxi, Smaxj;  // for local dynamic square memory
      int seq1begin, seq1end, seq2begin, seq2end;

      // alignment summary and detailed information
      // The following members are for holding secondary results
      int idencnt;
      int simcnt;
      int numgaps1;
      int numgaps2;
      int gaplen1;
      int gaplen2;

      //int alnlen;
      //vector<Alignedcol> alnresult;
      string topaln, middle, buttomaln;  // aligned sequences
      string topruler, buttomruler;
      ALNTYPE alntype;
};
class LSDynaln : public Dynaln {
   public:
      LSDynaln() : Dynaln() { }
      LSDynaln(const bioseq &s1, const bioseq &s2) 
         : Dynaln(s1,s2) { }
      LSDynaln(const bioseq &s1, const bioseq &s2, const string &mat) 
         : Dynaln(s1,s2,mat) { }
      ~LSDynaln();
      int global();
      int local();
      int runlocal() {
         int s=local(); buildResult(); return s; }
      //void global();
      /* globalLS in reverse direction
       */
      //int globalR(int b1, int e1, int b2, int e2);
      
      /* b1, e1, b2, and e2 are 0-based index of the 
       * first and the second sequences respectively.
       * if (b1<e1 and b2<e2)
       * then it does forward computation and store
       * the result in S.
       * when b1>e1 and b2>e2 then it does reverse
       * computation and store the result in SR
       *
       * return a pair of pointers to the top and buttom
       *    arrays respectively.
       *    This version is incorrect in detail.
       */
      //pair< int*, int*> computeScore(int b1, int e1, int b2, int e2);
      /* return the start index of the top and buttom rows
       * in the four matrices
       * This can be used to retrieve the top/buttom rows of the
       * three scoring matrices: M, IX, IY
       */
      pair<int, int> computeScoreFW(int b1, int e1, int b2, int e2);
      pair<int, int> computeScoreBW(int b1, int e1, int b2, int e2);

      int path(int b1, int e1, int b2, int e2);

      /* Completly different from the full dynamic program traceback.
       * Only the last two rows of the result matrix are in memory!
       *
       * given the two arrays M.first(TOP row), M.second (BUTTOM row)
       * and the 0-based index r, c in terms of the matrix
       * len1 is the length of the range of sequence 1
       * len2 is that of the sequence 2
       * if len1 and len2 are given, then it will use the 
       * reverse of the sequence for tracing.
       * Arguments: 
       *    r row index
       *    c column index, r and c specify the starting point
       *    for traceback and it is the cell with the best score.
       */
      //list<pair<int,int> > traceback(pair<int, int> mrow, int b1, int &e1, int b2, int &e2, bool rev=false);
      list<pair<int,int> > tracebackFW(pair<int,int> mrow, int b1, int &e1, int b2, int &e2);
      list<pair<int,int> > tracebackBW(pair<int,int> mrow, int b1, int &e1, int b2, int &e2);
      //void showaln();  // testing function

      void buildResult();

      void debug_showmatrixForward(pair<int,int> &rows, const int col);
      void debug_showmatrixBackward(pair<int,int> &rows, const int col);

   private:
      void allocmem();
      int *MR, *IXR, *IYR;
      //int *SR;  // for LOCAL this is used to memorize the traceback pointer
      //int *ItopR; // not used in local
};

#endif