alignment.cpp

   1 /*
   2  *  alignment.cpp
   3  *
   4  *  Created by Pat Schloss on 12/15/08.
   5  *  Copyright 2008 Patrick D. Schloss. All rights reserved.
   6  *
   7  *  This is a class for an abstract datatype for classes that implement various types of alignment      algorithms.
   8  *      As of 12/18/08 these included alignments based on blastn, needleman-wunsch, and the     Gotoh algorithms
   9  *
  10  */
  11
  12 #include "alignmentcell.hpp"
  13 #include "alignment.hpp"
  14
  15
  16 /**************************************************************************************************/
  17
  18 Alignment::Alignment() {        /*      do nothing      */      }
  19
  20 /**************************************************************************************************/
  21
  22 Alignment::Alignment(int A) : nCols(A), nRows(A) {
  23
  24         alignment.resize(nRows);                        //      For the Gotoh and Needleman-Wunsch we initialize the dynamic programming
  25         for(int i=0;i<nRows;i++){                       //      matrix by initializing a matrix that is A x A.  By default we will set A
  26                 alignment[i].resize(nCols);             //      at 2000 for 16S rRNA gene sequences
  27         }
  28
  29 }
  30
  31 /**************************************************************************************************/
  32
  33 void Alignment::traceBack(){                    //      This traceback routine is used by the dynamic programming algorithms
  34                                                                                 //      to fill the values of seqAaln and seqBaln
  35         seqAaln = "";
  36         seqBaln = "";
  37         int row = lB-1;
  38         int column = lA-1;
  39 //      seqAstart = 1;
  40 //      seqAend = column;
  41
  42         AlignmentCell currentCell = alignment[row][column];     //      Start the traceback from the bottom-right corner of the
  43                                                                                                                 //      matrix
  44
  45         if(currentCell.prevCell == 'x'){        seqAaln = seqBaln = "NOALIGNMENT";              }//If there's an 'x' in the bottom-
  46         else{                                                                                           //      right corner bail out because it means nothing got aligned
  47                 while(currentCell.prevCell != 'x'){                             //      while the previous cell isn't an 'x', keep going...
  48
  49                         if(currentCell.prevCell == 'u'){                        //      if the pointer to the previous cell is 'u', go up in the
  50                                 seqAaln = '-' + seqAaln;                                //      matrix.  this indicates that we need to insert a gap in
  51                                 seqBaln = seqB[row] + seqBaln;                  //      seqA and a base in seqB
  52                                 currentCell = alignment[--row][column];
  53                         }
  54                         else if(currentCell.prevCell == 'l'){           //      if the pointer to the previous cell is 'l', go to the left
  55                                 seqBaln = '-' + seqBaln;                                //      in the matrix.  this indicates that we need to insert a gap
  56                                 seqAaln = seqA[column] + seqAaln;               //      in seqB and a base in seqA
  57                                 currentCell = alignment[row][--column];
  58                         }
  59                         else{
  60                                 seqAaln = seqA[column] + seqAaln;               //      otherwise we need to go diagonally up and to the left,
  61                                 seqBaln = seqB[row] + seqBaln;                  //      here we add a base to both alignments
  62                                 currentCell = alignment[--row][--column];
  63                         }
  64                 }
  65         }
  66
  67         pairwiseLength = seqAaln.length();
  68         seqAstart = 1;  seqAend = 0;
  69         seqBstart = 1;  seqBend = 0;
  70
  71         for(int i=0;i<seqAaln.length();i++){
  72                 if(seqAaln[i] != '-' && seqBaln[i] == '-')              {       seqAstart++;    }
  73                 else if(seqAaln[i] == '-' && seqBaln[i] != '-') {       seqBstart++;    }
  74                 else                                                                                    {       break;                  }
  75         }
  76
  77         pairwiseLength -= (seqAstart + seqBstart - 2);
  78
  79         for(int i=seqAaln.length()-1; i>=0;i--){
  80                 if(seqAaln[i] != '-' && seqBaln[i] == '-')              {       seqAend++;              }
  81                 else if(seqAaln[i] == '-' && seqBaln[i] != '-') {       seqBend++;              }
  82                 else                                                                                    {       break;                  }
  83         }
  84         pairwiseLength -= (seqAend + seqBend);
  85
  86         seqAend = seqA.length() - seqAend - 1;
  87         seqBend = seqB.length() - seqBend - 1;
  88
  89 }
  90 /**************************************************************************************************/
  91
  92 Alignment::~Alignment(){
  93         try {
  94                 for (int i = 0; i < alignment.size(); i++) {
  95                         for (int j = (alignment[i].size()-1); j >= 0; j--) {  alignment[i].pop_back();  }
  96                 }
  97         }
  98         catch(exception& e) {
  99                 cout << "Standard Error: " << e.what() << " has occurred in the Alignment class Function ~Alignment. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 100                 exit(1);
 101         }
 102         catch(...) {
 103                 cout << "An unknown error has occurred in the Alignment class function ~Alignment. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 104                 exit(1);
 105         }
 106 }
 107
 108 /**************************************************************************************************/
 109
 110 string Alignment::getSeqAAln(){
 111         return seqAaln;                                                                         //      this is called to get the alignment of seqA
 112 }
 113
 114 /**************************************************************************************************/
 115
 116 string Alignment::getSeqBAln(){
 117         return seqBaln;                                                                         //      this is called to get the alignment of seqB
 118 }
 119
 120 /**************************************************************************************************/
 121
 122 int Alignment::getCandidateStartPos(){
 123         return seqAstart;                                                                       //      this is called to report the quality of the alignment
 124 }
 125
 126 /**************************************************************************************************/
 127
 128 int Alignment::getCandidateEndPos(){
 129         return seqAend;                                                                         //      this is called to report the quality of the alignment
 130 }
 131
 132 /**************************************************************************************************/
 133
 134 int Alignment::getTemplateStartPos(){
 135         return seqBstart;                                                                       //      this is called to report the quality of the alignment
 136 }
 137
 138 /**************************************************************************************************/
 139
 140 int Alignment::getTemplateEndPos(){
 141         return seqBend;                                                                         //      this is called to report the quality of the alignment
 142 }
 143
 144 /**************************************************************************************************/
 145
 146 int Alignment::getPairwiseLength(){
 147         return pairwiseLength;                                                          //      this is the pairwise alignment length
 148 }
 149
 150 /**************************************************************************************************/
 151
 152 //int Alignment::getLongestTemplateGap(){
 153 //
 154 //      int length = seqBaln.length();
 155 //      int longGap = 0;
 156 //      int gapLength = 0;
 157 //
 158 //      int start = seqAstart;
 159 //      if(seqAstart < seqBstart){      start = seqBstart;      }
 160 //      for(int i=seqAstart;i<length;i++){
 161 //              if(seqBaln[i] == '-'){
 162 //                      gapLength++;
 163 //              }
 164 //              else{
 165 //                      if(gapLength > 0){
 166 //                              if(gapLength > longGap){        longGap = gapLength;    }
 167 //                      }
 168 //                      gapLength = 0;
 169 //              }
 170 //      }
 171 //      return longGap;
 172 //}
 173
 174 /**************************************************************************************************/