overlap.cpp

   1 /*
   2  *  overlap.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 cleans up the alignment at the 3' end of the alignments.  Because the Gotoh and Needleman-Wunsch
   9  *      algorithms start the traceback from the lower-right corner of the dynamic programming matrix, there may be a lot of
  10  *      scattered bases in the alignment near the 3' end of the alignment.  Here we basically look for the largest score
  11  *      in the last column and row to determine whether there should be exta gaps in sequence A or sequence B.  The gap
  12  *      issues at the 5' end of the alignment seem to take care of themselves in the traceback.
  13  *
  14  */
  15 using namespace std;
  16
  17 #include "alignmentcell.hpp"
  18 #include "overlap.hpp"
  19
  20
  21 /**************************************************************************************************/
  22
  23 int Overlap::maxRow(vector<vector<AlignmentCell> >& alignment, const int band){
  24
  25         float max = -100;
  26         int end = lA - 1;
  27         int index = end;
  28
  29         for(int i=band;i<lB;i++){                                       //      find the row where the right most column has the highest alignment
  30                 if(alignment[i][end].cValue >= max){    //      score.
  31                         index = i;
  32                         max = alignment[i][end].cValue;
  33                 }
  34         }
  35         return index;
  36 }
  37
  38 /**************************************************************************************************/
  39
  40 int Overlap::maxColumn(vector<vector<AlignmentCell> >& alignment, const int band){
  41
  42         float max = -100;
  43         int end = lB - 1;
  44         int index = end;
  45
  46         for(int i=band;i<lA;i++){                                       //      find the column where the bottom most column has the highest
  47                 if(alignment[end][i].cValue >= max){    //      alignment score.
  48                         index = i;
  49                         max = alignment[end][i].cValue;
  50                 }
  51         }
  52         return index;
  53 }
  54
  55 /**************************************************************************************************/
  56
  57 void Overlap::setOverlap(vector<vector<AlignmentCell> >& alignment, const int nA, const int nB, const int band=0){
  58
  59         lA = nA;
  60         lB = nB;
  61
  62         int rowIndex = maxRow(alignment, band);         //      get the index for the row with the highest right hand side score
  63         int colIndex = maxColumn(alignment, band);      //      get the index for the column with the highest bottom row score
  64
  65         int row = lB-1;
  66         int column = lA-1;
  67
  68         if(colIndex == column && rowIndex == row){}     //      if the max values are the lower right corner, then we're good
  69         else if(alignment[row][colIndex].cValue < alignment[rowIndex][column].cValue){
  70                 for(int i=rowIndex+1;i<lB;i++){                 //      decide whether sequence A or B needs the gaps at the end either set
  71                         alignment[i][column].prevCell = 'u';//  the pointer upwards or...
  72                 }
  73
  74         }
  75         else {
  76                 for(int i=colIndex+1;i<lA;i++){
  77                         alignment[row][i].prevCell = 'l';       //      ...to the left
  78                 }
  79         }
  80 }                                                                                               //      the traceback should take care of the gaps at the 5' end
  81
  82 /**************************************************************************************************/
  83