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