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         try {
  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         catch(exception& e) {
  30                 errorOut(e, "Alignment", "Alignment");
  31                 exit(1);
  32         }
  33 }
  34 /**************************************************************************************************/
  35 void Alignment::resize(int A) {
  36         try {
  37                 nCols = A;
  38                 nRows = A;
  39
  40                 alignment.resize(nRows);
  41                 for(int i=0;i<nRows;i++){
  42                         alignment[i].resize(nCols);
  43                 }
  44         }
  45         catch(exception& e) {
  46                 errorOut(e, "Alignment", "resize");
  47                 exit(1);
  48         }
  49 }
  50 /**************************************************************************************************/
  51
  52 void Alignment::traceBack(){                    //      This traceback routine is used by the dynamic programming algorithms
  53         try {
  54                                                                         //      to fill the values of seqAaln and seqBaln
  55                 seqAaln = "";
  56                 seqBaln = "";
  57                 int row = lB-1;
  58                 int column = lA-1;
  59                 //      seqAstart = 1;
  60                 //      seqAend = column;
  61
  62                 AlignmentCell currentCell = alignment[row][column];     //      Start the traceback from the bottom-right corner of the
  63                 //      matrix
  64
  65                 if(currentCell.prevCell == 'x'){        seqAaln = seqBaln = "NOALIGNMENT";              }//If there's an 'x' in the bottom-
  66                 else{                                                                                           //      right corner bail out because it means nothing got aligned
  67                         while(currentCell.prevCell != 'x'){                             //      while the previous cell isn't an 'x', keep going...
  68
  69                                 if(currentCell.prevCell == 'u'){                        //      if the pointer to the previous cell is 'u', go up in the
  70                                         seqAaln = '-' + seqAaln;                                //      matrix.  this indicates that we need to insert a gap in
  71                                         seqBaln = seqB[row] + seqBaln;                  //      seqA and a base in seqB
  72                                         currentCell = alignment[--row][column];
  73                                 }
  74                                 else if(currentCell.prevCell == 'l'){           //      if the pointer to the previous cell is 'l', go to the left
  75                                         seqBaln = '-' + seqBaln;                                //      in the matrix.  this indicates that we need to insert a gap
  76                                         seqAaln = seqA[column] + seqAaln;               //      in seqB and a base in seqA
  77                                         currentCell = alignment[row][--column];
  78                                 }
  79                                 else{
  80                                         seqAaln = seqA[column] + seqAaln;               //      otherwise we need to go diagonally up and to the left,
  81                                         seqBaln = seqB[row] + seqBaln;                  //      here we add a base to both alignments
  82                                         currentCell = alignment[--row][--column];
  83                                 }
  84                         }
  85                 }
  86
  87                 pairwiseLength = seqAaln.length();
  88                 seqAstart = 1;  seqAend = 0;
  89                 seqBstart = 1;  seqBend = 0;
  90
  91                 for(int i=0;i<seqAaln.length();i++){
  92                         if(seqAaln[i] != '-' && seqBaln[i] == '-')              {       seqAstart++;    }
  93                         else if(seqAaln[i] == '-' && seqBaln[i] != '-') {       seqBstart++;    }
  94                         else                                                                                    {       break;                  }
  95                 }
  96
  97                 pairwiseLength -= (seqAstart + seqBstart - 2);
  98
  99                 for(int i=seqAaln.length()-1; i>=0;i--){
 100                         if(seqAaln[i] != '-' && seqBaln[i] == '-')              {       seqAend++;              }
 101                         else if(seqAaln[i] == '-' && seqBaln[i] != '-') {       seqBend++;              }
 102                         else                                                                                    {       break;                  }
 103                 }
 104                 pairwiseLength -= (seqAend + seqBend);
 105
 106                 seqAend = seqA.length() - seqAend - 1;
 107                 seqBend = seqB.length() - seqBend - 1;
 108         }
 109         catch(exception& e) {
 110                 errorOut(e, "Alignment", "traceBack");
 111                 exit(1);
 112         }
 113 }
 114 /**************************************************************************************************/
 115
 116 Alignment::~Alignment(){
 117         try {
 118                 for (int i = 0; i < alignment.size(); i++) {
 119                         for (int j = (alignment[i].size()-1); j >= 0; j--) {  alignment[i].pop_back();  }
 120                 }
 121         }
 122         catch(exception& e) {
 123                 errorOut(e, "Alignment", "~Alignment");
 124                 exit(1);
 125         }
 126 }
 127
 128 /**************************************************************************************************/
 129
 130 string Alignment::getSeqAAln(){
 131         return seqAaln;                                                                         //      this is called to get the alignment of seqA
 132 }
 133
 134 /**************************************************************************************************/
 135
 136 string Alignment::getSeqBAln(){
 137         return seqBaln;                                                                         //      this is called to get the alignment of seqB
 138 }
 139
 140 /**************************************************************************************************/
 141
 142 int Alignment::getCandidateStartPos(){
 143         return seqAstart;                                                                       //      this is called to report the quality of the alignment
 144 }
 145
 146 /**************************************************************************************************/
 147
 148 int Alignment::getCandidateEndPos(){
 149         return seqAend;                                                                         //      this is called to report the quality of the alignment
 150 }
 151
 152 /**************************************************************************************************/
 153
 154 int Alignment::getTemplateStartPos(){
 155         return seqBstart;                                                                       //      this is called to report the quality of the alignment
 156 }
 157
 158 /**************************************************************************************************/
 159
 160 int Alignment::getTemplateEndPos(){
 161         return seqBend;                                                                         //      this is called to report the quality of the alignment
 162 }
 163
 164 /**************************************************************************************************/
 165
 166 int Alignment::getPairwiseLength(){
 167         return pairwiseLength;                                                          //      this is the pairwise alignment length
 168 }
 169
 170 /**************************************************************************************************/
 171
 172 //int Alignment::getLongestTemplateGap(){
 173 //
 174 //      int length = seqBaln.length();
 175 //      int longGap = 0;
 176 //      int gapLength = 0;
 177 //
 178 //      int start = seqAstart;
 179 //      if(seqAstart < seqBstart){      start = seqBstart;      }
 180 //      for(int i=seqAstart;i<length;i++){
 181 //              if(seqBaln[i] == '-'){
 182 //                      gapLength++;
 183 //              }
 184 //              else{
 185 //                      if(gapLength > 0){
 186 //                              if(gapLength > longGap){        longGap = gapLength;    }
 187 //                      }
 188 //                      gapLength = 0;
 189 //              }
 190 //      }
 191 //      return longGap;
 192 //}
 193
 194 /**************************************************************************************************/