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() { }
      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()
            { 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()
            { 
               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();

      /* This methods build addional information from the
       * initial phase of the algorithm.
       *    1. call traceback.
       *    2. clear existing result (call clearResult())
       *    3. then use the alnidx list to update other information:
       *       idencnt: counts for identical resudues
       *       simcnt: similar residues
       *       numgaps1, numgaps2: number of gaps in seq1 and seq2
       *       gaplen1, gaplen2: total gaps in both sequences.
       *       topaln: aligned string of seq1
       *       buttomaln: aligned string of seq2
       *       middle: alignment middle line (| for identical : for similar)
       * This method create additional redundant information for
       * easy usage and user interface.
       */
      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; }
      string toDelimitedString(const string &dl=",") 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(); }
      /* count the number of gaps in the two sequences
       */
      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;
      /* use & | operator to figure out which direction
       * store the traceback pointers in an integer
       * multiple pointers can be stored as sum of the 
       * individual pointers
       * DIAG+TOP=3 etc...
       */
      static const int PTRNULL=0;
      static const int PTRDIAG=1;
      static const int PTRTOP=2; // Ix
      static const int PTRLEFT=4; // Iy

      /* debug functions */
      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);
      /* this version works with the pointers
       * and only affect the alnidx list.
       */
      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;
      
      /* alnidx contains final alignment result
       * the numbers are the index in the two sequences.
       * -1 for a gap
       * we can print this result in many different ways
       */
      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
      /* the two numbers are use to record the
       * traceback starting point in terms of index
       * in the input sequences.
       * For global, it is seq1len-1, seq2len-1
       */
      int Smaxi, Smaxj;  // for local dynamic square memory

      // alignment summary and detailed information
      int idencnt;
      int simcnt;
      int numgaps1;
      int numgaps2;
      int gaplen1;
      int gaplen2;

      /* start and end of alignment, if local then
       * it is the aligned part, if global, it 
       * excludes the begining and ending gaps
       * 0-based index, inclusive
      */
      int seq1begin, seq1end, seq2begin, seq2end;
      //int alnlen;
      //vector<Alignedcol> alnresult;
      string topaln, middle, buttomaln;  // aligned sequences
      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();
      /* the linear space version
       * b1: begin index of sequence1 (0-based index)
       * e1: one-passed the end index of sequence1 (0-based index)
       * return: the best score
       */
      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);
      /* works on the reverse matrices
       * b1>e1 | b2>e2
       */
      pair<int, int> computeScoreBW(int b1, int e1, int b2, int e2);

      /* arguments:
         * 0-based index of sequences
         *  b1, e1:  begin and end of sequence 1
         *  b2, e2: begin and end of sequence 2
         *  Here b1<e1 and b2<e2
         *
         *  The result will be consolidated in the alnidx list.
       */
      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);
      /* arguments: 
       *    mrow  index of the last two rows of the matrix (top,buttom)
       *    This function will decrement e1 and e2.  If they have become
       *    smaller than b1 and b2 respectively, then the job has been
       *    finished and there is no need to do another recursion call.
       *    The caller (path()) function uses this information.
       * return: path in forward direction
       */
      list<pair<int,int> > tracebackFW(pair<int,int> mrow, int b1, int &e1, int b2, int &e2);
      /*
       * arguments:
       *    mrow is the matrix row, a pair of integers for the top and buttom
       *       array starting point.  This information was used to locate the
       *       starting point for the MR, IXR, and IYR arrays.  The first
       *       integer starts the top array, and the second integer starts the
       *       buttom array.  MR + mrow.first --> top array.
       * Return the path in the forward direction, there is no need to reverse it
       */
      list<pair<int,int> > tracebackBW(pair<int,int> mrow, int b1, int &e1, int b2, int &e2);
      //void showaln();  // testing function

      /* totally overwrite the base class version.
       * There are some commonality, need to be factored out
       * in the future.
       */
      void buildResult();

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

   private:
      /* allocate memory use two rows for top problem
       * two rows for buttom problem
       * S and SR have only one row to store the 
       * score of the last row.
       */
      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