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
  16 #include "alignmentcell.hpp"
  17 #include "overlap.hpp"
  18
  19
  20 /**************************************************************************************************/
  21
  22 int Overlap::maxRow(vector<vector<AlignmentCell> >& alignment, const int band){
  23
  24         float max = -100;
  25         int end = lA - 1;
  26         int index = end;
  27
  28         for(int i=band;i<lB;i++){                                       //      find the row where the right most column has the highest alignment
  29                 if(alignment[i][end].cValue >= max){    //      score.
  30                         index = i;
  31                         max = alignment[i][end].cValue;
  32                 }
  33         }
  34         return index;
  35 }
  36
  37 /**************************************************************************************************/
  38
  39 int Overlap::maxColumn(vector<vector<AlignmentCell> >& alignment, const int band){
  40
  41         float max = -100;
  42         int end = lB - 1;
  43         int index = end;
  44
  45         for(int i=band;i<lA;i++){                                       //      find the column where the bottom most column has the highest
  46                 if(alignment[end][i].cValue >= max){    //      alignment score.
  47                         index = i;
  48                         max = alignment[end][i].cValue;
  49                 }
  50         }
  51         return index;
  52 }
  53
  54 /**************************************************************************************************/
  55
  56 void Overlap::setOverlap(vector<vector<AlignmentCell> >& alignment, const int nA, const int nB, const int band=0){
  57
  58         lA = nA;
  59         lB = nB;
  60
  61         int rowIndex = maxRow(alignment, band);         //      get the index for the row with the highest right hand side score
  62         int colIndex = maxColumn(alignment, band);      //      get the index for the column with the highest bottom row score
  63
  64         int row = lB-1;
  65         int column = lA-1;
  66
  67         if(colIndex == column && rowIndex == row){}     //      if the max values are the lower right corner, then we're good
  68         else if(alignment[row][colIndex].cValue < alignment[rowIndex][column].cValue){
  69                 for(int i=rowIndex+1;i<lB;i++){                 //      decide whether sequence A or B needs the gaps at the end either set
  70                         alignment[i][column].prevCell = 'u';//  the pointer upwards or...
  71                 }
  72
  73         }
  74         else {
  75                 for(int i=colIndex+1;i<lA;i++){
  76                         alignment[row][i].prevCell = 'l';       //      ...to the left
  77                 }
  78         }
  79 }                                                                                               //      the traceback should take care of the gaps at the 5' end
  80
  81 /**************************************************************************************************/
  82