QueryManager/MatchAligner.h

Go to the documentation of this file.

/*  Last edited: May 26 15:54 2002 (ac2) */

// #######################################################################

// SSAHA : Sequence Search and Alignment by Hashing Algorithm
// Version 3.2, released 1st March 2004
// Copyright (c) Genome Research 2002

// SSAHA is free software; you can redistribute it and/or modify 
// it under the terms of version 2 of the GNU General Public Licence
// as published by the Free Software Foundation.
 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public Licence for more details.
 
// You should have received a copy of the GNU General Public Licence
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
// or see the on-line version at http://www.gnu.org/copyleft/gpl.txt

// #######################################################################

// Module Name  : MatchAligner
// File Name    : MatchAligner.h
// Language     : C++
// Module Author: Anthony J. Cox (ac2@sanger.ac.uk)

// Include guard:
#ifndef INCLUDED_MatchAligner
#define INCLUDED_MatchAligner

// Description:

// Includes:
#include "MatchStore.h"
#include "SequenceEncoder.h"
#include <cassert>

// NB it is good practise for #include statements in header files to be
// replaced by forward declarations if at all possible

// ### Class Declarations ###

// Declarations for MatchTaskAlign

typedef unsigned char uchar;

typedef unsigned char AlignBreakType;

const AlignBreakType alignBreakNull = 0;
const AlignBreakType alignBreakMismatch = 1;
const AlignBreakType alignBreakGapQuery = 2;
const AlignBreakType alignBreakGapSubject = 3;
const AlignBreakType alignBreakFrameShiftQuery = 4;
const AlignBreakType alignBreakFrameShiftSubject = 5;

typedef int ScoreType;
typedef ScoreType ScoreTable[65536];

const ScoreType matchScoreDNA    = 1; // | matchMask;
const ScoreType mismatchScoreDNA = -2;
const ScoreType gapStartScoreDNA = -4;
const ScoreType gapExtendScoreDNA  = -3;
const ScoreType nScoreDNA = 0; // same as crossmatch defaults
const ScoreType veryBadScoreIndeed = -1000000000;


// next 3 are as used by NCBI - see URL in .h file
const ScoreType stopCodonScoreBLOSUM = -4; 
const ScoreType gapStartScoreBLOSUM = -16; // = -31; // = -10;
const ScoreType gapExtendScoreBLOSUM =-1;

// frameShiftScoreBLOSUM is set so that three frame shifts incur less penalty
// than shifting by a codon
// Algorithm should then be able to pick up alignments of the form
// A..-B..-C..-D.. <-- three frame shifts, score -15
// ||| ||| ||| |||
// nnnnnnnnnnnnnnn 
// rather than just doing
// A..B..C..---D.. <-- these count as a single gap, score -16
// |||xxxxxx   |||
// nnnnnnnnnnnnnnn
// Also frame shift score ought to be worse than worse mismatch score
const ScoreType frameShiftScoreBLOSUM = -5; // = -10; // -3;


// from new alignment alg
//const ScoreType veryBadScoreIndeed = -1000000000;
//const ScoreType gapScore = -4;
//const ScoreType matchScore = 1;
//const ScoreType mismatchScore = -2; // -10
// const ScoreType frameShiftScore = -3;
typedef unsigned char PathType;
const PathType fromN = (PathType)0x1;
const PathType fromW = (PathType)0x2;
const PathType fromSW = (PathType)0x4;
const PathType fromPrevFrame1 = (PathType)0x8;
const PathType fromPrevFrame2 = (PathType)0x10;
const PathType fromFinished = (PathType)0x20;  

// ----------------------------------------------
// Declarations for MatchTaskAlign and subclasses
// ----------------------------------------------

class MatchAligner;

class MatchTaskAlign : public MatchTask
{
public:
  MatchTaskAlign( SourceReader& querySource, 
                  SourceReader& subjectSource,
                  MatchAligner* pAlign,
                  bool reverseQueryCoords=true, 
                  bool doAlignment=true,
                  ostream& outputStream=cout );
  virtual ~MatchTaskAlign();

  virtual void operator()(MatchStore& store);

  void reverseQueryCoords( bool yesOrNo= true) 
  { 
    reverseQueryCoords_ = yesOrNo;
  } // ~void reverseQueryCoords( bool yesOrNo= true) 

  void doAlignment( bool yesOrNo= true) 
  { 
    doAlignment_ = yesOrNo;
  } // ~void doAlignment( bool yesOrNo= true) 
 

protected:

  MatchAligner* pAlign_;
  ostream& outputStream_;
  SourceReader& querySource_;
  SourceReader& subjectSource_;
  // reverseQueryCoords_ - if `true' reverses orientation of
  // query coordinates, so they are as seen in the forward frame
  // Basically set to true for ssaha executable but false in 
  // ssahaClient, because in that case any necessary reversals
  // are done in ssahaServer
  bool reverseQueryCoords_;
  // doAlignment_ - if `true' prints match description line plus
  // graphical alignment for each match, if `false' just prints
  // match description line
  bool doAlignment_;

}; // ~class MatchTaskAlign : public MatchTask



// --------------------------------------------
// Declarations for MatchAligner and subclasses
// --------------------------------------------

// alignment of two sequences is described by a vector of AlignInfo 
// Meant for sequences that contain large stretches of matches
// punctuated by the occasional AlignBreak

struct AlignInfo
{
  AlignInfo( void ) : numMatches(0), breakType(0) {}
  int numMatches;
  AlignBreakType breakType;
}; // ~struct AlignInfo

struct Alignment : public deque<AlignInfo>
{
  ScoreType totalScore_;
}; // ~struct Alignment



// Class Name : StaticScoreTable
// Description: A static instance of this is created for each
// score table required. Implements `singleton' design pattern.
class StaticScoreTable
{
  typedef void (* makeTablePointer)( ScoreTable& table );
 public:
  StaticScoreTable( makeTablePointer p ) : makeTable_(p) {}
  ScoreTable* getTable() // TBD const ptr 
  { 
    if (makeTable_) (*makeTable_)(data_); 
    makeTable_ = NULL; 
    return &data_;
  }
 private:
  ScoreTable data_;
  makeTablePointer makeTable_;
}; // ~class StaticScoreTable


void makeScoreTableDNA( ScoreTable& table );
void makeScoreTableBlosum62( ScoreTable& table );





class MatchAligner
{
public:
  MatchAligner(   int numCols,
                  int bandExtension,
                  ScoreTable* pTable,
                  ostream& outputStream );
  virtual ~MatchAligner();

  void operator()( const char* pQuery, int queryStart, int queryEnd,
                   const char* pSubject, int subjectStart, int subjectEnd );
  virtual void createAlignment
  ( Alignment& alignment,
    const char* pQuery, int queryStart, int queryEnd,
    const char* pSubject, int subjectStart, int subjectEnd );

  virtual void formatAlignment
  ( Alignment& alignment,
    const char* pQuery, int queryStart, int queryEnd,
    const char* pSubject, int subjectStart, int subjectEnd );

  bool matchChar( const char a, const char b )
  {
    return (tolower(a)==tolower(b));
  }

  void outputAlignmentColumn
    ( const char queryChar, const char AlignChar, const char subjectChar,
      int queryNext, int subjectNext);

  void outputAlignmentLine( void );

protected:

  ostream& outputStream_;

  // These store incomplete lines of formatted alignment while
  // awaiting output
  char* pBufSeq1_;
  char* pBufAlign_;
  char* pBufSeq2_;

  // These point to the current characters in the line buffers above
  char* pCursorSeq1_;
  char* pCursorAlign_;
  char* pCursorSeq2_;

  // This is the number of characters per formatted line
  int numCols_;
  int bandExtension_;

  static StaticScoreTable tableDNA_;
  static StaticScoreTable tableBlosum62_;

  ScoreTable* pTable_;


}; // ~class MatchAligner


// Class Name: MatchAlignerDNA
// Description: Aligns using DNA score table
class MatchAlignerDNA : public MatchAligner
{
 public:
  MatchAlignerDNA(  int numCols=80,
                    int bandExtension=0,
                    ostream& outputStream=cout ) :
  MatchAligner( numCols, 
                bandExtension,
                tableDNA_.getTable(),
                outputStream )
    {} // ~MatchTaskAlignDNA::MatchTaskAlignDNA
}; // ~class MatchAlignerDNA : public MatchAligner


// Class Name: MatchAlignerProtein
// Description: Aligns protein using score table based on BLOSUM62
class MatchAlignerProtein : public MatchAligner
{
 public:
  MatchAlignerProtein( int numCols=80,
                       int bandExtension=0,
                       ostream& outputStream=cout ) :
    MatchAligner( numCols, 
                  bandExtension,
                  tableBlosum62_.getTable(), 
                  outputStream )
    {} // ~MatchAlignerProtein::MatchAlignerProtein

}; // ~class MatchAlignerProtein : public MatchAligner

// Class Name: MatchAlignerTranslated
// Description: 
class MatchAlignerTranslated : public MatchAligner
{
 public: // TBD protected
  MatchAlignerTranslated(   int numCols,
                            int bandExtension,
                            bool isQueryProtein,
                            bool isSubjectProtein,
                            ostream& outputStream );

  void codonize( const char* pSeq, 
                 int seqSize,
                 int finalFrame,
                 vector<vector<char> >& translatedSeqs ); 

  virtual void createAlignment
  ( Alignment& alignment,
    const char* pQuery, int queryStart, int queryEnd,
    const char* pSubject, int subjectStart, int subjectEnd );

  static char getCodon( const char* pChar )
  {
    return (  (    (ttDNA[ *pChar ] ==nv)
                   || (ttDNA[ *(pChar+1) ] ==nv)
                   || (ttDNA[ *(pChar+2) ] ==nv) )
              ? 'X'
              : gResidueNames[ ttCodon[ ttDNA[ *(pChar) ] << 4
                                      | ttDNA[ *(pChar+1) ] << 2
                                      | ttDNA[ *(pChar+2) ] ] ] );
  } // ~getCodon


 protected:
  bool isQueryProtein_;
  bool isSubjectProtein_;
};


// Class Name: MatchAlignerTranslatedProtein
// Description: Aligns protein and a DNA sequence.
// If isQueryProtein assumes query is protein, subject is DNA, 
// else vice versa.
class MatchAlignerTranslatedProtein : public MatchAlignerTranslated
{
 public:
  MatchAlignerTranslatedProtein
    ( int isQueryProtein = true,
      int numCols=80,
      int bandExtension=0,
      ostream& outputStream=cout );

  virtual void formatAlignment
  ( Alignment& alignment,
    const char* pQuery, int queryStart, int queryEnd,
    const char* pSubject, int subjectStart, int subjectEnd );

 protected:
}; // ~class MatchAlignerTranslatedProtein : public MatchAlignerTranslated

// Class Name: MatchAlignerTranslatedDNA
// Description: Aligns translations of two DNA sequences. 
class MatchAlignerTranslatedDNA : public MatchAlignerTranslated
{
 public:
  MatchAlignerTranslatedDNA( int numCols=80,
                             int bandExtension=0,
                             ostream& outputStream=cout ) :
    MatchAlignerTranslated( numCols, 
                            bandExtension,
                            false, 
                            false, 
                            outputStream )
    {} // ~MatchAlignerTranslatedDNA::MatchAlignerTranslatedDNA

  virtual void formatAlignment
  ( Alignment& alignment,
    const char* pQuery, int queryStart, int queryEnd,
    const char* pSubject, int subjectStart, int subjectEnd );
}; // ~class MatchAlignerTranslatedDNA : public MatchAlignerProteinDNA

// --------------------------------------------
// Data for scoring matches between amino acids
// --------------------------------------------

#ifdef DATA_SOURCES_FOR_BLOSUM62_MATRIX
The below is cut and pasted from
http://www.ncbi.nlm.nih.gov/IEB/ToolBox/C_DOC/lxr/source/data/BLOSUM62

-- quote --

  1 #  Matrix made by matblas from blosum62.iij
  2 #  * column uses minimum score
  3 #  BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
  4 #  Blocks Database = /data/blocks_5.0/blocks.dat
  5 #  Cluster Percentage: >= 62
  6 #  Entropy =   0.6979, Expected =  -0.5209
  7    A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  *
  8 A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1  0 -4 
  9 R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1  0 -1 -4 
 10 N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  3  0 -1 -4 
 11 D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4  1 -1 -4 
 12 C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4 
 13 Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0  3 -1 -4 
 14 E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4 
 15 G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -2 -1 -4 
 16 H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0  0 -1 -4 
 17 I -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3 -3 -3 -1 -4 
 18 L -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1 -4 -3 -1 -4 
 19 K -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2  0  1 -1 -4 
 20 M -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1 -3 -1 -1 -4 
 21 F -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1 -3 -3 -1 -4 
 22 P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2 -2 -1 -2 -4 
 23 S  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2  0  0  0 -4 
 24 T  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0 -1 -1  0 -4 
 25 W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3 -4 -3 -2 -4 
 26 Y -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1 -3 -2 -1 -4 
 27 V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4 -3 -2 -1 -4 
 28 B -2 -1  3  4 -3  0  1 -1  0 -3 -4  0 -3 -3 -2  0 -1 -4 -3 -3  4  1 -1 -4 
 29 Z -1  0  0  1 -3  3  4 -2  0 -3 -3  1 -1 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4 
 30 X  0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2  0  0 -2 -1 -1 -1 -1 -1 -4 
 31 * -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4  1 

-- unquote --

In addition:
Score of any character with - is -10
Score of - with - is -1
These are the gap open and gap extension penalties used with BLOSUM62 in
Huang/Zhang, CABIOS 12(6) 1996, pp. 497-506

#endif

static const char blosumResidues[] 
= "ARNDCQEGHILKMFPSTWYVBZX";
// 12345678901234567890123

static const ScoreType blosumScores[] =
{
4, // A
-1, 5, // R
-2, 0, 6, // N
-2, -2, 1, 6, // D
0, -3, -3, -3, 9,// C
-1, 1, 0, 0, -3, 5,// Q
-1, 0, 0, 2, -4, 2, 5,// E
0, -2, 0, -1, -3, -2, -2, 6, // G
-2, 0, 1, -1, -3, 0, 0, -2, 8, // H
-1, -3, -3, -3, -1, -3, -3, -4, -3, 4, // I
-1, -2, -3, -4, -1, -2, -3, -4, -3, 2, 4, // L
-1, 2, 0, -1, -3, 1, 1, -2, -1, -3, -2, 5, // K
-1, -1, -2, -3, -1, 0, -2, -3, -2, 1, 2, -1, 5, // M
-2, -3, -3, -3, -2, -3, -3, -3, -1, 0, 0, -3, 0, 6, // F
-1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4, 7, // P
1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2, 0, -1, -2, -1, 4, // S
0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 1, 5, // T
-3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1, 1, -4, -3, -2, 11, // W
-2, -2, -2, -3, -2, -1, -2, -3, 2, -1, -1, -2, -1, 3, -3, -2, -2, 2, 7,// Y
0, -3, -3, -3, -1, -2, -2, -3, -3, 3, 1, -2, 1, -1, -2, -2, 0, -3, -1, 4,// V
-2, -1, 3, 4, -3, 0, 1, -1, 0, -3, -4, 0, -3, -3, -2, 0, -1, -4, -3, -3, 4,// B
-1, 0, 0, 1, -3, 3, 4, -2, 0, -3, -3, 1, -1, -3, -1, 0, -1, -3, -2, -2, 1, 4, // Z
0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, 0, 0, -2, -1, -1, -1, -1, -1 //X
};

static const int numBlosumResidues = 23;


// ---------------------------------------------------------
// Declarations for standard (2D) banded dynamic programming
// ---------------------------------------------------------


// Class name : PathMatrix
// Description: Idea of PathMatrix is that all 4 types of banded
// dynamic programming needed can be handled by using different
// types for PATH_TYPE
// 2D (`standard') banded d.p.: DNA vs DNA, protein vs protein
// 3D banded d.p.: protein vs translated DNA
// 4D banded d.p.: translated DNA vs translated DNA
// COLUMN_FILLER fills in cells column by column
// TRACE_BACKER computes the traceback path cell by cell

// If the standard dp matrix is ...
//
//     -   p1  p1  p1  p1  p1 
//  -  c00 c01 c02 c03 c04 c05
// q1  c10 c11 c12 c13 c14 c15
// q2  c20 c21 c22 c23 c24 c25
// q3  c30 c31 c32 c33 c34 c35
// q4  c40 c41 c42 c43 c44 c45
// q5  c50 c51 c52 c53 c54 c55
// q6  c60 c61 c62 c63 c64 c65
// q7  c70 c71 c72 c73 c74 c75
// q8  c80 c81 c82 c83 c84 c85
//
// ... then the corresponding banded matrix with band extension 0 is ...
//
// c00 c11 c22 c33 c44 c55
// c10 c21 c32 c43 c54 c65
// c20 c31 c42 c53 c64 c75
// c30 c41 c52 c63 c74 c85
//
//  ... and the corresponding banded matrix with band extension 2 is ...
//
// xxx xxx c02 c12 c24 c35
// xxx c01 c12 c23 c34 c45
// c00 c11 c22 c33 c44 c55
// c10 c21 c32 c43 c54 c65
// c20 c31 c42 c53 c64 c75
// c30 c41 c52 c63 c74 c85
// c40 c51 c62 c73 c84 xxx
// c50 c61 c72 c83 xxx xxx

// Notes:
// 1. Each column is a vector<PathType> of size 
// length(q1)-length(p1)+1+(2*bandExtension)
// This means if bandExtension is nonzero the elements marked xxx
// are never used, even though storage is allocated.
// 2. When doing the traceback, cell traversals W NW N in the standard matrix
// correspond to cell traversals SW W N respectively in the banded matrix.
// 3. Only directional pointers are stored in each cell. Only the current
// and previous column of scores are retained as the cells are filled in.

template <class PATH_TYPE> class PathMatrix 
  : public vector<vector<PATH_TYPE> >
{
public:
  typedef pair<vector<vector<PATH_TYPE> >::iterator,
    vector<PATH_TYPE>::iterator> CellIterator;

  template<class MATRIX_FILLER> ScoreType fillIn( MATRIX_FILLER& doMatrix )
  {
    return doMatrix(*this);
  } // ~fillIn

  template <class TRACE_BACKER> void traceBack
  ( TRACE_BACKER& doCell )
  {
    // NB we are tracing back from finish to start
    //    CellIterator start; 
    //    start.first = &front();
    //  start.second = front().begin()+startPos;
    //    finish.first =  &back();  
    //    finish.second = back().begin()+finishPos;
 
    CellIterator current = lastCell_;
    //    current.first =  &back();  
    //  current.second = back().begin()+finishPos;
    //    cout << "Cell iterator - column: " << current.first-begin()
    //  << " row: " << current.second-current.first->begin() << endl;
    while( doCell(current) );
    // {
    //   cout << "Cell iterator - column: " << current.first-begin()
    //   << " row: " << current.second-current.first->begin() << endl;
    // }

  } // ~traceBack

  CellIterator lastCell_;

}; // ~class PathMatrix

// Class name : ScoreMaker
// Description: Returns score of two chars according to scoreTable
class ScoreMaker
{
 public:
  ScoreMaker( const ScoreTable& scoreTable ) :
    scoreTable_(scoreTable),
    b_((unsigned short*)&a_[0])
    {}
  ScoreType operator()( uchar c1, uchar c2 )
  {
    a_[0]=c1;
    a_[1]=c2;
    return scoreTable_[ *b_ ];
  }
  ScoreType getGapScore( void )
  {
    a_[0]='a'; // a is both a valid nucleotide and amino acid code
    a_[1]='-';
    return scoreTable_[ *b_ ];
  }
 private:
  const ScoreTable& scoreTable_;
  uchar a_[2];
  unsigned short* b_;
}; // ~class ScoreMaker


// Class name : CellFiller
// Description:
class CellFiller
{
public:
  CellFiller( void ) :  max_(0) {}
  ScoreType operator()
  ( PathType& cell, 
    ScoreType scoreW,
    ScoreType scoreSW,
    ScoreType scoreN )
    {
      max_=max( max( scoreW, scoreSW), scoreN );
      cell |= fromW * (scoreW==max_);
      cell |= fromSW * (scoreSW==max_);
      cell |= fromN * (scoreN==max_);
      //      cout << "Cell: " 
      //   << ((cell & fromW)?'-':'.')
      //   << ((cell & fromSW)?'/':'.') 
      //   << ((cell & fromN)?'|':'.') 
      //   << " " << scoreW
      //   << " " << scoreSW
      //   << " " << scoreN 
      //   << " " << max_ << endl; 
      return max_;
    } // ~CellFiller::operator()

private:
  ScoreType max_;
};

// Class name : ColumnFillerBasic
// Description: This produces a column of a PathMatrix 
// for DNA vs DNA or protein vs protein matching
class ColumnFillerBasic
{
public:
  ColumnFillerBasic( const char* p1, int p1Size, 
                     const char* p2, int p2Size,
                     int bandExtension,
                     const ScoreTable& scoreTable) :
    p1_(p1), 
    p2_(p2), 
    //    bandExtension_(bandExtension),
    bandExtension_( (bandExtension<(p1Size/2))
                    ? bandExtension
                    : max( (p1Size/2)-1, 0 ) ), 
    bandWidth_(p2Size-p1Size+1),
    bandLength_(p1Size+1),
    colSize_(p2Size-p1Size+1+(2*bandExtension_)),
    //    scoresPossible_(3),
    fillCell_(),
    v1_(colSize_, veryBadScoreIndeed ), 
    v2_(colSize_, veryBadScoreIndeed ),
    pLast_(&v1_),
    pCurrent_(&v2_),
    getScore_(scoreTable),
    gapScore_(getScore_.getGapScore())
    {
      // TBD check bandExtension_ is not too big
    }
  ScoreType operator()( PathMatrix<PathType>& matrix );
private:
  const char* p1_;
  const char* p2_;
  int bandExtension_;
  const int bandWidth_;
  const int bandLength_;
  const int colSize_;
  CellFiller fillCell_;
  //  vector<pair<ScoreType, PathType> > scoresPossible_;
  vector<ScoreType> v1_;
  vector<ScoreType> v2_;
  vector<ScoreType>* pLast_;
  vector<ScoreType>* pCurrent_;
  vector<ScoreType>* temp_;
  ScoreMaker getScore_;
  ScoreType  gapScore_;

}; // ~class ColumnFillerBasic

// Class name : TraceBackerBasic
// Description: This carries out one step of the traceback
// for DNA vs DNA or protein vs protein matching
class TraceBackerBasic
{
public:
  TraceBackerBasic( AlignBreakType p1Gap, AlignBreakType p2Gap,
                    Alignment& alignment ) :
    p1Gap_(p1Gap), p2Gap_(p2Gap),
    alignment_(alignment) {}

  bool operator()( PathMatrix<PathType>::CellIterator& current );
 private:
  AlignBreakType p1Gap_;
  AlignBreakType p2Gap_;
  Alignment& alignment_;
}; // ~class TraceBackerBasic


// ------------------------------------------------------------
// Definitions for banded dynamic programming with translations
// ------------------------------------------------------------

template< typename T> class Array2D
{
public:
  Array2D() { data_[0]=0; data_[1]=0; data_[2]=0; }
  Array2D( T i, T j, T k) 
    { data_[0]=i; data_[1]=j; data_[2] = k; }
  T& operator[]( unsigned int i ) { return data_[i]; }
private:
  T data_[gNumReadingFrames];
};

template< typename T> class Array3D
{
public:
  Array3D() {}
  Array3D( const Array2D<T>& i, 
           const Array2D<T>& j, 
           const Array2D<T>& k) 
    { data_[0]=i; data_[1]=j; data_[2] = k; }
  Array2D<T>& operator[]( unsigned int i ) { return data_[i]; }
private:
  Array2D<T> data_[gNumReadingFrames];
};

const Array2D<ScoreType> veryBadScore2D( veryBadScoreIndeed, 
                              veryBadScoreIndeed, 
                              veryBadScoreIndeed );

const Array3D<ScoreType> veryBadScore3D( veryBadScore2D, 
                                         veryBadScore2D, 
                                         veryBadScore2D );


typedef Array3D<ScoreType> ScoreType3D;
typedef Array3D<PathType> PathType3D;


// Class name : CellFiller3D
// Description:
class CellFiller3D
{
public:
  CellFiller3D( void ) :  max_(0) {}
  ScoreType operator()
  ( PathType& cell, 
    ScoreType scoreW,
    ScoreType scoreSW,
    ScoreType scoreN,
    ScoreType scoreFrame1,
    ScoreType scoreFrame2 )
  {

    max_=max( max( max( scoreW, 
                        scoreSW ), 
                   max( scoreFrame1, 
                        scoreFrame2 ) ),
              scoreN );
    cell |= fromW * (scoreW==max_);
    cell |= fromSW * (scoreSW==max_);
    cell |= fromN * (scoreN==max_);
    cell |= fromPrevFrame1 * (scoreFrame1==max_);
    cell |= fromPrevFrame2 * (scoreFrame2==max_);
    //    cout << "Cell: " 
    //   << ((cell & fromW)?'-':'.')
    //     << ((cell & fromSW)?'/':'.') 
    //   << ((cell & fromN)?'|':'.') 
    //    << ((cell & fromPrevFrame1)?'<':'.') 
    //   << ((cell & fromPrevFrame2)?'^':'.') 
    //   << " " << scoreW
    //   << " " << scoreSW
    //   << " " << scoreN 
    //   << " " << scoreFrame1
    //   << " " << scoreFrame2
    //   << " " << max_ << endl; 
    return max_;
  } // ~CellFiller3D::operator()


private:
  ScoreType max_;
};


// Class name : ColumnFiller3D
// Description: This produces a column of a PathMatrix 
// for protein vs translated DNA or translated DNA vs translated DNA matching
class ColumnFiller3D
{
public:
  //  ColumnFiller3D( const char* p1, int p1Size, bool isProtein1, 
  //      const char* p2, int p2Size, bool isProtein2,
  //      int bandExtension,
  //      const ScoreTable& scoreTable );
  ColumnFiller3D( const char** p1Trans, int p1Size, int p1FinalFrame,
                  const char** p2Trans, int p2Size, int p2FinalFrame,
                  int bandExtension,
                  const ScoreTable& scoreTable );

  ScoreType operator()( PathMatrix<PathType3D>& matrix );

  // NB all these template functions are specialized (ie overridden)
  // for the case false. This allows them to be `switched off' at
  // compile time. Definitions are in MatchAligner.cpp

  template< bool canGoW> ScoreType getScoreW
  ( int j, int k, int l )
  {
    return (*pLast_)[j][k][l];
  } // ~template< bool canGoW> ScoreType getScoreW

  template< bool canGoSW> ScoreType getScoreSW
  ( int j, int k, int l )
  {
    return (*pLast_)[j+1][k][l];
  } // ~template< bool canGoSW> ScoreType getScoreSW

  template< bool canGoN> ScoreType getScoreN
  ( int j, int k, int l )
  {
    return (*pCurrent_)[j-1][k][l];
  } // ~template< bool canGoN> ScoreType getScoreN


  // Function: getScorePrevFrame1
  // If the frame k is 1 or 2, the score for the previous frame
  // is in the same cell. 
  // If the frame k is 0, need to grab the frame 2 score of the
  // cell to the SW. This corresponds to the W-ward cell in the
  // full d.p. matrix
  template< bool canGoSW > ScoreType getScorePrevFrame1
  ( int j, int k, int l )
  {
    return ((k!=0) 
            ? (*pCurrent_)[j][k-1][l] 
            : ( (numFrames1_==1)
                ? veryBadScoreIndeed
                : getScoreSW<canGoSW>( j, 2, l) ) );
  } // ~template< bool canGoSW > ScoreType getScorePrevFrame1

  // Function: getScorePrevFrame2
  // If the frame l is 1 or 2, the score for the previous frame
  // is in the same cell. 
  // If the frame l is 0, need to grab the frame 2 score of the
  // cell to the N. This corresponds to the N-wards cell in the
  // full d.p. matrix
  template< bool canGoN > ScoreType getScorePrevFrame2
  ( int j, int k, int l )
  {
    return ((l!=0) 
            ? (*pCurrent_)[j][k][l-1] 
            : ( (numFrames2_==1)
                ? veryBadScoreIndeed
                :getScoreN<canGoN>( j, k, 2) ) );
  } // ~template< bool canGoN > ScoreType getScorePrevFrame2

  template< bool canMatch> ScoreType getScoreChar
  ( int i, int j, int k, int l)
  {
    assert(p1_[k]!=NULL);
    assert(p2_[l]!=NULL);

    //        cout << "getScoreChar "
    //  << p1_[k][i-1] << p2_[l][j-bandExtension_+i-1] << " "
    //  << i << " " 
    //  << j << " " 
    //  << k << " " 
    //  << l << endl;
    return (getScore_(p1_[k][i-1],p2_[l][j-bandExtension_+i-1]));
  } // ~template< bool canMatch> ScoreType getScoreChar



  template< bool canGoW, bool canGoSW, bool canGoN, bool canMatch> 
  void fillCell3D
  ( PathMatrix<PathType3D>& matrix, int i, int j )
  {
    //        cout << "filling cell " << i << " " << j << " : "
    //  << (canGoW?'T':'F')
    //  << (canGoSW?'T':'F')
    //  << (canGoN?'T':'F')
    //   << (canMatch?'T':'F')
    //  << endl;

    for ( int k(0); k < numFrames1_; k++ )
    {
      for ( int l(0); l < numFrames2_; l++ )
      {
        //              cout << k << " " << l << endl;
        (*pCurrent_)[j][k][l] 
          = fillCell_( matrix[i][j][k][l],
                       getScoreW<canGoW>( j, k, l)
                       + getScoreChar<canMatch>(i,j,k,l),
                       getScoreSW<canGoSW>( j, k, l)
                       + gapStartScoreBLOSUM,
                       getScoreN<canGoN>( j, k, l)
                       + gapStartScoreBLOSUM,
                       getScorePrevFrame1<canGoSW>( j, k, l)
                       + frameShiftScoreBLOSUM,
                       getScorePrevFrame2<canGoN>( j, k, l) 
                       + frameShiftScoreBLOSUM );
      } // ~for l
    } // ~for k
  } // ~fillCell3D



private:
  const char* p1_[gNumReadingFrames];
  const char* p2_[gNumReadingFrames];
  vector< vector<char> > p1Translations_;
  vector< vector<char> > p2Translations_;

  int bandExtension_;
  int bandWidth_;
  int bandLength_;
  int colSize_;
  int finalFrame1_;
  int finalFrame2_;
  int numFrames1_;
  int numFrames2_;

  CellFiller3D fillCell_;
  ScoreMaker getScore_;
  //  vector<pair<ScoreType, PathType> > scoresPossible_;
  vector<ScoreType3D> v1_;
  vector<ScoreType3D> v2_;
  vector<ScoreType3D>* pLast_;
  vector<ScoreType3D>* pCurrent_;
  vector<ScoreType3D>* temp_;



}; // ~class ColumnFiller3D


// Class name : TraceBacker3D
// Description: This carries out one step of the traceback
// for matches involving translations or protein
class TraceBacker3D
{
public:
  //  TraceBacker3D( AlignBreakType p1Gap, AlignBreakType p2Gap,
  //                int p1FinalFrame, int p2FinalFrame,
  //                Alignment& alignment ) :
  //  p1Gap_(p1Gap), p2Gap_(p2Gap),
  //  k_(p1FinalFrame), 
  //  l_(p2FinalFrame),
  //  alignment_(alignment) {}
  TraceBacker3D( int queryFinalFrame, int subjectFinalFrame,
                 bool p1IsQuery,
                 Alignment& alignment ) :
    p1Gap_( (p1IsQuery) ? 
            alignBreakGapQuery : alignBreakGapSubject ),
    p2Gap_( (p1IsQuery) ? 
            alignBreakGapSubject : alignBreakGapQuery ),
    p1FrameShift_( (p1IsQuery) ? 
                   alignBreakFrameShiftQuery : alignBreakFrameShiftSubject ),
    p2FrameShift_( (p1IsQuery) ? 
                   alignBreakFrameShiftSubject : alignBreakFrameShiftQuery ),
    k_( (p1IsQuery) ? 
        queryFinalFrame : subjectFinalFrame ), 
    l_( (p1IsQuery) ? 
        subjectFinalFrame : queryFinalFrame ), 
    alignment_(alignment) 
    {}
  bool operator()( PathMatrix<PathType3D>::CellIterator& current );
private:

  AlignBreakType p1Gap_;
  AlignBreakType p2Gap_;
  AlignBreakType p1FrameShift_;
  AlignBreakType p2FrameShift_;

  int k_; // frame of p1
  int l_; // frame of p2
  Alignment& alignment_;



};




// ### Function Declarations ###

// Name:
// Arguments:
// TYPE  NAME  IN/OUT COMMENT
// Returns: TYPE COMMENT

// End of include guard:
#endif

// End of file MatchAligner.h

---