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