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