needlemanoverlap.cpp

   1 /*
   2  *  needleman.cpp
   3  *
   4  *
   5  *  Created by Pat Schloss on 12/15/08.
   6  *  Copyright 2008 Patrick D. Schloss. All rights reserved.
   7  *
   8  *      This class is an Alignment child class that implements the Gotoh pairwise alignment algorithm as described in:
   9  *
  10  *              Gotoh O. 1982.  An improved algorithm for matching biological sequences.  J. Mol. Biol.  162:705-8.
  11  *              Myers, EW & Miller, W.  1988.  Optimal alignments in linear space.  Comput Appl Biosci. 4:11-7.
  12  *
  13  *      This method is nice because it allows for an affine gap penalty to be assessed, which is analogous to what is used
  14  *      in blast and is an alternative to Needleman-Wunsch, which only charges the same penalty for each gap position.
  15  *      Because this method typically has problems at the ends when two sequences do not full overlap, we employ a separate
  16  *      method to fix the ends (see Overlap class documentation)
  17  *
  18  */
  19
  20 using namespace std;
  21
  22 #include "alignmentcell.hpp"
  23 #include "alignment.hpp"
  24 #include "overlap.hpp"
  25 #include "needlemanoverlap.hpp"
  26
  27 /**************************************************************************************************/
  28
  29 NeedlemanOverlap::NeedlemanOverlap(float gO, float m, float mm, int r) ://      note that we don't have a gap extend
  30 gap(gO), match(m), mismatch(mm), Alignment(r) {                                                 //      the gap openning penalty is assessed for
  31                                                                                                                                                 //      every gapped position
  32         for(int i=1;i<nCols;i++){
  33                 alignment[0][i].prevCell = 'l';                                 //      initialize first row by pointing all poiters to the left
  34                 alignment[0][i].cValue = 0;                                             //      and the score to zero
  35         }
  36
  37         for(int i=1;i<nRows;i++){
  38                 alignment[i][0].prevCell = 'u';                                 //      initialize first column by pointing all poiters upwards
  39                 alignment[i][0].cValue = 0;                                             //      and the score to zero
  40         }
  41 }
  42
  43 /**************************************************************************************************/
  44
  45 NeedlemanOverlap::~NeedlemanOverlap(){  /*      do nothing      */      }
  46
  47 /**************************************************************************************************/
  48
  49 void NeedlemanOverlap::align(string A, string B){
  50
  51         seqA = ' ' + A; lA = seqA.length();             //      algorithm requires a dummy space at the beginning of each string
  52         seqB = ' ' + B; lB = seqB.length();             //      algorithm requires a dummy space at the beginning of each string
  53
  54         for(int i=1;i<lB;i++){                                  //      This code was largely translated from Perl code provided in Ex 3.1
  55                 for(int j=1;j<lA;j++){                          //      of the O'Reilly BLAST book.  I found that the example output had a
  56                                                                                         //      number of errors
  57                         float diagonal;
  58                         if(seqB[i] == seqA[j])  {       diagonal = alignment[i-1][j-1].cValue + match;          }
  59                         else                                    {       diagonal = alignment[i-1][j-1].cValue + mismatch;       }
  60
  61                         float up        = alignment[i-1][j].cValue + gap;
  62                         float left      = alignment[i][j-1].cValue + gap;
  63
  64                         if(diagonal >= up){
  65                                 if(diagonal >= left){
  66                                         alignment[i][j].cValue = diagonal;
  67                                         alignment[i][j].prevCell = 'd';
  68                                 }
  69                                 else{
  70                                         alignment[i][j].cValue = left;
  71                                         alignment[i][j].prevCell = 'l';
  72                                 }
  73                         }
  74                         else{
  75                                 if(up >= left){
  76                                         alignment[i][j].cValue = up;
  77                                         alignment[i][j].prevCell = 'u';
  78                                 }
  79                                 else{
  80                                         alignment[i][j].cValue = left;
  81                                         alignment[i][j].prevCell = 'l';
  82                                 }
  83                         }
  84                 }
  85         }
  86         Overlap over;
  87         over.setOverlap(alignment, lA, lB, 0);          //      Fix gaps at the beginning and end of the sequences
  88         traceBack();                                                            //      Traceback the alignment to populate seqAaln and seqBaln
  89 }
  90
  91 /**************************************************************************************************/
  92