QueryManager/MatchAligner.cpp

Go to the documentation of this file.

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

// 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.cpp
// Language     : C++
// Module Author: Anthony J. Cox (ac2@sanger.ac.uk)

// Description:

// Includes:

#include "MatchAligner.h"
#include "SequenceReader.h"
#include "SequenceEncoder.h"
#include <iomanip> // for setprecision() etc...

// ### Function Definitions ###

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

// Definitions for MatchTaskAlign

// MatchTaskAlign function definitions

MatchTaskAlign::MatchTaskAlign
( SourceReader& querySource, 
  SourceReader& subjectSource,
  MatchAligner* pAlign,
  bool reverseQueryCoords,
  bool doAlignment,
  ostream& outputStream ):
  querySource_(querySource), 
  subjectSource_(subjectSource),
  pAlign_(pAlign),
  //  numCols_(numCols),
  reverseQueryCoords_(reverseQueryCoords),
  doAlignment_(doAlignment),
  outputStream_(outputStream)
{
} // ~MatchTaskAlign::MatchTaskAlign


MatchTaskAlign::~MatchTaskAlign()
{
  delete pAlign_;
} // ~MatchTaskAlign::~MatchTaskAlign


void MatchTaskAlign::operator()(MatchStore& store)
{
  //  cout << "MTA::()" << endl;

  if (store.empty()) return;

  SequenceOffset effectiveQueryStart, effectiveQueryEnd;


  vector<Match*>::const_iterator i(store.begin());

  //  outputStream_ << endl << "Matches For Query " 
  //       << (*i)->getQueryNum()
  //       << " (" << (*i)->getQuerySize()
  //       << " bases): " << (*i)->getQueryName()
  //       << "\n\n";    

  outputStream_ << setprecision(2) << setiosflags(ios::fixed);

  //  vector<char> queryData, subjectData;

  for( ;i!=store.end();++i)
  {
    
    if( reverseQueryCoords_ && (!(*i)->isQueryForward()) )
    {
      effectiveQueryStart = (*i)->getQuerySize()-(*i)->getQueryEnd()+1;
      effectiveQueryEnd = (*i)->getQuerySize()-(*i)->getQueryStart()+1;
    } // ~if
    else
    {
      effectiveQueryStart = (*i)->getQueryStart();
      effectiveQueryEnd = (*i)->getQueryEnd();
    } // ~else

    outputStream_ 
        << ((*i)->isQueryForward()?"F":"R") 
        << ((*i)->isSubjectForward()?"F":"R") << "\t"
        << (*i)->getQueryName() << "\t"
        << effectiveQueryStart << "\t"     
        << effectiveQueryEnd << "\t"     
        << (*i)->getSubjectName() << "\t"
        << (*i)->getSubjectStart() << "\t"
        << (*i)->getSubjectEnd() << "\t"
        << (*i)->getNumBases() << "\t" 
        << 100.0*((*i)->getNumBases()) / 
      ((*i)->getQueryEnd()-(*i)->getQueryStart()+1)
        << endl;

    // suppress output of graphical alignments if required
    if( !doAlignment_ ) continue;


    char* pSubject;
    char* pQuery;
  
    //    cout << "extracting subject" << endl;
    
    if ((*i)->isSubjectForward())
    {
      //      cout << "extracting subject fwd" << endl;
      subjectSource_.extractSource( &pSubject, //subjectData, 
                                    (*i)->getSubjectNum(),
                                    (*i)->getSubjectStart(),
                                    (*i)->getSubjectEnd() );
    } // ~if
    else
    {
      //      cout << "extracting subject rev" << endl;
      subjectSource_.extractSourceReverse( &pSubject, //subjectData, 
                                    (*i)->getSubjectNum(),
                                    (*i)->getSubjectStart(),
                                    (*i)->getSubjectEnd() );
    } // ~else

    if (pSubject[0]=='\0') continue; // did not get subject data from server

    if ((*i)->isQueryForward())
    {
      //      cout << "extracting query fwd" << endl;
      querySource_.extractSource
        ( &pQuery, //queryData, 
          (*i)->getQueryNum(),
          effectiveQueryStart,
          effectiveQueryEnd );
    } // ~if
    else
    {
      //      cout << "extracting query rev" << endl;
      querySource_.extractSourceReverse
      ( &pQuery, //queryData, 
        (*i)->getQueryNum(),
        effectiveQueryStart,  
        effectiveQueryEnd );
    } // ~else
    //   cout << "all extractions done" << endl;

    (*pAlign_)( pQuery,//(const char*) &queryData[0], 
           effectiveQueryStart,  
           effectiveQueryEnd,
           pSubject,//(const char*) &subjectData[0],
           (*i)->getSubjectStart(),
           (*i)->getSubjectEnd() );

    //    queryData.clear();
    //    subjectData.clear();
  } // ~for i


} // ~MatchTaskAlign::operator()(MatchStore& store)


// MatchAligner function definitions


StaticScoreTable MatchAligner::tableDNA_(&makeScoreTableDNA);
StaticScoreTable MatchAligner::tableBlosum62_(&makeScoreTableBlosum62);


MatchAligner::MatchAligner
( int numCols,
  int bandExtension,
  ScoreTable* pTable,
  ostream& outputStream ) :
  numCols_(numCols),
  bandExtension_(bandExtension),
  pTable_(pTable),
  outputStream_(outputStream)
{
  pBufSeq1_= new char [numCols+1];
  pBufSeq2_= new char [numCols+1];
  pBufAlign_= new char [numCols+1];
  pBufSeq1_[numCols]='\0';
  pBufSeq2_[numCols]='\0';
  pBufAlign_[numCols]='\0';
} // ~MatchAligner::MatchAligner


MatchAligner::~MatchAligner()
{
  delete [] pBufSeq1_;
  delete [] pBufSeq2_;
  delete [] pBufAlign_;
} // ~MatchAligner::~MatchAligner



void MatchAligner::operator()
( const char* pQuery, int queryStart, int queryEnd,
  const char* pSubject, int subjectStart, int subjectEnd )
{

  Alignment alignment;
  createAlignment
    ( alignment, 
      pQuery, queryStart, queryEnd, 
      pSubject, subjectStart, subjectEnd );

  formatAlignment
    ( alignment, 
      pQuery, queryStart, queryEnd, 
      pSubject, subjectStart, subjectEnd );

} // ~void MatchTaskAlign::align


void MatchAligner::createAlignment
( 
 Alignment& alignment,
 const char* pQuery, int queryStart, int queryEnd,
 const char* pSubject, int subjectStart, int subjectEnd )
{

  const char* p1;
  const char* p2;

  int querySize(queryEnd-queryStart+1);
  int subjectSize(subjectEnd-subjectStart+1);

  //  int bandWidth;
  // int bandLength;
  int p1Size, p2Size;

  //  AlignBreakType gapInLargest, gapInSmallest;
  AlignBreakType p1Gap, p2Gap;

  // ensure *p1 <= *p2
  if (querySize<=subjectSize)
  {
    p1=pQuery; p2=pSubject;
    p1Size=querySize; p2Size = subjectSize;
    //    bandWidth  = subjectSize-querySize+1;
    //   bandLength = querySize;
    p1Gap=alignBreakGapQuery;
    p2Gap=alignBreakGapSubject;
  } // ~if
  else
  {
    p1=pSubject; p2=pQuery;
    p1Size= subjectSize; p2Size = querySize;
    //    bandWidth  = querySize-subjectSize+1;
    //   bandLength = subjectSize;
    p1Gap=alignBreakGapSubject;
    p2Gap=alignBreakGapQuery;
  } // ~else

  PathMatrix<PathType> path;
  ColumnFillerBasic doMatrix( p1, p1Size, p2, p2Size, bandExtension_, 
                              *pTable_ ); 
  // 0 = band ext
  alignment.totalScore_ = path.fillIn(doMatrix);
  // print(path);
  TraceBackerBasic doCell(p1Gap, p2Gap, alignment);
  path.traceBack(doCell); // 0 = band ext
  //  path.traceBack(doCell, bandExtension_, 
  //     p2Size-p1Size+bandExtension_); // 0 = band ext


} // ~void MatchAligner::createAlignment


void MatchAligner::formatAlignment
( 
 Alignment& alignment,
 const char* pQuery, int queryStart, int queryEnd,
 const char* pSubject, int subjectStart, int subjectEnd )
{

  outputStream_ << "Alignment score: " << alignment.totalScore_ << endl;

  int queryNext(queryStart);
  int subjectNext(subjectStart);

  pCursorSeq1_ = &pBufSeq1_[numCols_]; // triggers reset to start of line

  *pBufSeq1_='\0'; // ensures 3 null strings are printed at first line break 
  //  *pBufSeq2_='\0'; // - this line not necessary??
  //  *pBufAlign_='\0'; // - this line not necessary??

  Alignment::iterator i(alignment.begin());

  while( i!=alignment.end() ) 
  {
    //    cout << queryNext << " " << subjectNext << endl;
    
    if (i->numMatches>0)
    { // TBD do these en masse?
      //      outputAlignmentColumn
      //        (*(pQuery++), '|', *(pSubject++), queryNext++, subjectNext++ ); 
      outputAlignmentColumn
      ( *pQuery, 
        ((tolower(*pQuery)==tolower(*pSubject))?'|':'x'), 
        *pSubject, 
        queryNext++, 
        subjectNext++ ); 
      pQuery++;
      pSubject++;

      //   queryNext++;
      //  subjectNext++;
    i->numMatches--;
    }
    else
    {
      if (i->breakType==alignBreakMismatch)
      {
        outputAlignmentColumn
          (*(pQuery++), 'x', *(pSubject++), queryNext++, subjectNext++ ); 
        //      queryNext++;
        //      subjectNext++;
      } // ~if
      else if (i->breakType==alignBreakGapQuery)
      {
        outputAlignmentColumn
          ( '-', ' ', *(pSubject++), queryNext, subjectNext++ ); 
        //      subjectNext++;
      } // ~else if
      else if (i->breakType==alignBreakGapSubject)
      {
        outputAlignmentColumn
          (*(pQuery++), ' ', '-', queryNext++, subjectNext ); 
        //      queryNext++;
      } // ~else if
      i++;
    } // ~else

  } // ~while

  *pCursorSeq1_='\0';
  *pCursorSeq2_='\0';
  *pCursorAlign_='\0';

  outputAlignmentLine();

} // ~void MatchAligner::formatAlignment

void MatchAligner::outputAlignmentColumn
( const char queryChar, const char alignChar, const char subjectChar,
  int queryNext, int subjectNext)
{
    if (pCursorSeq1_==&pBufSeq1_[numCols_])
    {
      if (*pBufSeq1_!='\0')
        outputAlignmentLine();
      sprintf( pBufSeq1_, "Q:%9.9d ", queryNext );
      sprintf( pBufAlign_, "            " ); // should be 12 spaces
      sprintf( pBufSeq2_, "S:%9.9d ", subjectNext );
      pCursorSeq1_ = &pBufSeq1_[12];
      pCursorSeq2_ = &pBufSeq2_[12];
      pCursorAlign_ = &pBufAlign_[12];

    } // ~if
  *(pCursorSeq1_++)  = queryChar;
  *(pCursorAlign_++) = alignChar;
  *(pCursorSeq2_++)  = subjectChar;
} // ~void MatchAligner::outputAlignmentColumn




void MatchAligner::outputAlignmentLine( void )
{
  outputStream_ << pBufSeq1_ << endl << pBufAlign_ << endl 
                << pBufSeq2_ << endl << endl;
} // ~void MatchTaskAlign::outputAlignmentLine( void )



void makeScoreTableDNA( ScoreTable& table )
{

  uchar a[2];
  unsigned short* b = (unsigned short*) &a[0];
  //  a[0]=0; a[1]=0;

  const char base[] = {"agct"};

  ScoreType thisScore;
  for (int i(0); i < 4; i++ )
  {
    for (int j(0); j < 4; j++ )
    {
      thisScore = ((i==j) ? matchScoreDNA: mismatchScoreDNA);
      a[0]=base[i]; a[1] = base[j]; table[ *b ] = thisScore;
      a[1] = toupper(base[j]); table[ *b ] = thisScore;
      a[0]=toupper(base[i]); a[1] = base[j]; table[ *b ] = thisScore;
      a[1] = toupper(base[j]); table[ *b ] = thisScore;
    } // ~for j

    a[0]=base[i];  
    a[1] = 'n'; table[ *b ] = nScoreDNA;
    a[1] = 'N'; table[ *b ] = nScoreDNA;
    a[1] = '-'; table[ *b ] = gapStartScoreDNA;

    a[0]=toupper(base[i]);  
    a[1] = 'n'; table[ *b ] = nScoreDNA;
    a[1] = 'N'; table[ *b ] = nScoreDNA;
    a[1] = '-'; table[ *b ] = gapStartScoreDNA;

    a[1]=base[i];  
    a[0] = 'n'; table[ *b ] = nScoreDNA;
    a[0] = 'N'; table[ *b ] = nScoreDNA;
    a[0] = '-'; table[ *b ] = gapStartScoreDNA;

    a[1]=toupper(base[i]);  
    a[0] = 'n'; table[ *b ] = nScoreDNA;
    a[0] = 'N'; table[ *b ] = nScoreDNA;
    a[0] = '-'; table[ *b ] = gapStartScoreDNA;
    


  } // ~for i;

    a[0]='n'; a[1] = 'n'; table[ *b ] = nScoreDNA;
    a[0]='n'; a[1] = 'N'; table[ *b ] = nScoreDNA;
    a[0]='N'; a[1] = 'n'; table[ *b ] = nScoreDNA;
    a[0]='N'; a[1] = 'N'; table[ *b ] = nScoreDNA;

    a[0]='-'; a[1] = 'n'; table[ *b ] = gapStartScoreDNA;
    a[0]='-'; a[1] = 'N'; table[ *b ] = gapStartScoreDNA;
    a[0]='n'; a[1] = '-'; table[ *b ] = gapStartScoreDNA;
    a[0]='N'; a[1] = '-'; table[ *b ] = gapStartScoreDNA;

    a[0]='-'; a[1] = '-'; table[ *b ] = gapExtendScoreDNA;

} // ~void makeScoreTableDNA( ScoreTable& table )


void makeScoreTableBlosum62( ScoreTable& table )
{

  uchar a[2];
  unsigned short* b = (unsigned short*) &a[0];

  const ScoreType* pScore = blosumScores;


  for (int i(0); i < numBlosumResidues; i++ )
  {
    for (int j(0); j < i; j++ )
    {
      //      cout << i << " " << j << endl; 

      // a[0] = uc, a[1] = uc
      a[0]=blosumResidues[i]; a[1]=blosumResidues[j]; table[ *b ] = *pScore;
      //      cout << a[0] << a[1] << " " table[ *b ] << endl;
      // a[0] = uc, a[1] = lc
      a[1]=tolower(a[1]); table[ *b ] = *pScore;
      // a[0] = lc, a[1] = lc
      a[0]=tolower(a[0]); table[ *b ] = *pScore;
      // a[0] = lc, a[1] = uc
      a[1]=toupper(a[1]); table[ *b ] = *pScore;

      a[0]=blosumResidues[j]; a[1]=blosumResidues[i]; table[ *b ] = *pScore;
      a[1]=tolower(a[1]); table[ *b ] = *pScore;
      a[0]=tolower(a[0]); table[ *b ] = *pScore;
      a[1]=toupper(a[1]); table[ *b ] = *pScore;
     
      ++pScore;
    } // ~for j

    a[0] = blosumResidues[i];
    a[1] = blosumResidues[i]; table[ *b ] = (*pScore); // a0=uc, a1=uc
    a[1] = tolower(a[1]); table[ *b ] = (*pScore);  // a0=uc, a1=lc
    a[1] = '*'; table[ *b ] = stopCodonScoreBLOSUM; // a0=uc, a1=*
    a[1] = '-'; table[ *b ] = gapStartScoreBLOSUM;  // a0=uc, a1=-

    a[0] = tolower(a[0]); 
    a[1] = blosumResidues[i]; table[ *b ] = (*pScore); // a0=lc, a1=uc
    a[1] = tolower(a[1]); table[ *b ] = (*pScore);  // a0=lc, a1=lc
    a[1] = '*'; table[ *b ] = stopCodonScoreBLOSUM; // a0=lc, a1=*
    a[1] = '-'; table[ *b ] = gapStartScoreBLOSUM;  // a0=lc, a1=-

    a[1] = blosumResidues[i];
    a[0] = '*'; table[ *b ] = stopCodonScoreBLOSUM; // a0=*, a1=uc
    a[0] = '-'; table[ *b ] = gapStartScoreBLOSUM;  // a0=-, a1=uc

    a[1] = tolower(a[1]);
    a[0] = '*'; table[ *b ] = stopCodonScoreBLOSUM; // a0=*, a1=lc
    a[0] = '-'; table[ *b ] = gapStartScoreBLOSUM;  // a0=-, a1=lc

    ++pScore;
  } // ~for i

  a[0] = '*'; a[1] = '*'; table[ *b ] = 1;
  a[0] = '*'; a[1] = '-'; table[ *b ] = gapStartScoreBLOSUM;
  a[0] = '-'; a[1] = '*'; table[ *b ] = gapStartScoreBLOSUM;
  a[0] = '-'; a[1] = '-'; table[ *b ] = gapExtendScoreBLOSUM;

} // ~void makeScoreTableBlosum62( ScoreTable& table )

MatchAlignerTranslated::MatchAlignerTranslated
(   int numCols,
    int bandExtension,
    bool isQueryProtein,
    bool isSubjectProtein,
    ostream& outputStream ) :
  MatchAligner( numCols, 
                bandExtension,
                tableBlosum62_.getTable(), 
                outputStream ),
  isQueryProtein_(isQueryProtein),
  isSubjectProtein_(isSubjectProtein)
{

} // ~MatchAlignerTranslated::MatchAlignerTranslated


void MatchAlignerTranslated::codonize
( const char* pSeq, 
  int seqSize,
  int finalFrame,
  vector<vector<char> >& translatedSeqs )
{
  translatedSeqs.resize(gNumReadingFrames);
  translatedSeqs[0].resize(seqSize);
  translatedSeqs[1].resize(seqSize - (finalFrame==0));
  translatedSeqs[2].resize(seqSize - (finalFrame!=2));
  //  cout << "codonizing: " << seqSize << endl;
  const char* pChar;
  vector<char>::iterator i;

  for (int j(0); j < gNumReadingFrames; j++ )
  {
    pChar = pSeq+j;
    for (i=translatedSeqs[j].begin(); 
         i!=translatedSeqs[j].end(); i++, 
         pChar+=gNumReadingFrames )
    {    
      //       cout << *pChar << *(pChar+1) << *(pChar+2) << endl;
      //  cout << int ( ttDNA[ *(pChar++ ] << 4
      //     | (ttDNA[ *(pChar++) ] << 2)
      //     | ttDNA[ *(pChar++) ] ) << endl; 
      if (    (ttDNA[ *(pChar)   ]==nv)
           || (ttDNA[ *(pChar+1) ]==nv)
           || (ttDNA[ *(pChar+2) ]==nv) )
      { 
        *i='X';
      } // ~if   
      else
      {
        *i= gResidueNames[ ttCodon[ ttDNA[ *(pChar) ] << 4
                                  | ttDNA[ *(pChar+1) ] << 2
                                  | ttDNA[ *(pChar+2) ] ] ];
      } // ~else

    } // ~for i

  } // ~for j

} // ~MatchAlignerTranslated::codonize


void MatchAlignerTranslated::createAlignment
( Alignment& alignment,
  const char* pQuery, int queryStart, int queryEnd,
  const char* pSubject, int subjectStart, int subjectEnd )
{
  PathMatrix<PathType3D> path;

  int queryFinalFrame(0);
  int subjectFinalFrame(0);

  int querySize(queryEnd-queryStart+1);
  int subjectSize(subjectEnd-subjectStart+1);

  vector<vector<char> > queryTranslated;
  vector<vector<char> > subjectTranslated;
  const char* pQueryTrans[gNumReadingFrames] = 
  { pQuery, NULL, NULL };
  const char* pSubjectTrans[gNumReadingFrames] =
  { pSubject, NULL, NULL };

  //  AlignBreakType p1Gap, p2Gap;
  //  int p1FinalFrame, p2FinalFrame;
  bool p1IsQuery;


  if (!isQueryProtein_)
  {
    queryFinalFrame = querySize % gNumReadingFrames;
    querySize /= gNumReadingFrames;
    codonize( pQuery, querySize, queryFinalFrame, queryTranslated );
    queryTranslated[1].push_back(0);
    queryTranslated[2].push_back(0);
    pQueryTrans[0] = static_cast<const char*>(&*queryTranslated[0].begin());
    pQueryTrans[1] = static_cast<const char*>(&*queryTranslated[1].begin());
    pQueryTrans[2] = static_cast<const char*>(&*queryTranslated[2].begin());
  } // ~if

  if (!isSubjectProtein_)
  {
    subjectFinalFrame = subjectSize % gNumReadingFrames;
    subjectSize /= gNumReadingFrames;
    codonize( pSubject, subjectSize, subjectFinalFrame, subjectTranslated );
    subjectTranslated[1].push_back(0);
    subjectTranslated[2].push_back(0);
    pSubjectTrans[0]= static_cast<const char*>(&*subjectTranslated[0].begin());
    pSubjectTrans[1]= static_cast<const char*>(&*subjectTranslated[1].begin());
    pSubjectTrans[2]= static_cast<const char*>(&*subjectTranslated[2].begin());
  } // ~if

  if (querySize<=subjectSize)
  { // then p1 = query, p2 = subject
    p1IsQuery = true;
    ColumnFiller3D doMatrix( pQueryTrans, querySize, queryFinalFrame,
                             pSubjectTrans, subjectSize, subjectFinalFrame,
                             bandExtension_, *pTable_ ); 
    alignment.totalScore_ = path.fillIn(doMatrix);

  } // ~if
  else
  { // p1 = subject, p2 = query
    p1IsQuery = false;
    ColumnFiller3D doMatrix( pSubjectTrans, subjectSize, subjectFinalFrame,
                             pQueryTrans, querySize, queryFinalFrame,
                             bandExtension_, *pTable_ ); 
    alignment.totalScore_ = path.fillIn(doMatrix);

  } // ~else

  TraceBacker3D doCell( queryFinalFrame, subjectFinalFrame, 
                        p1IsQuery, alignment );
  path.traceBack(doCell); 
  //  path.traceBack(doCell, bandExtension_, 
  //     fabs(subjectSize-querySize)+bandExtension_); // 0 = band ext

  if (!isQueryProtein_)
  {
    queryTranslated[1].pop_back();
    queryTranslated[2].pop_back();
  }
  if (!isSubjectProtein_)
  {
    subjectTranslated[1].pop_back();
    subjectTranslated[2].pop_back();
  }

} // ~void MatchAlignerTranslated::createAlignment


MatchAlignerTranslatedProtein::MatchAlignerTranslatedProtein
( int isQueryProtein,
  int numCols,
  int bandExtension,
  ostream& outputStream ) :
  MatchAlignerTranslated( numCols, 
                          bandExtension,
                          isQueryProtein, 
                          !isQueryProtein, 
                          outputStream )
{
} // ~MatchAlignerTranslatedProtein::MatchAlignerTranslatedProtein




void MatchAlignerTranslatedProtein::formatAlignment
( Alignment& alignment,
  const char* pQuery, int queryStart, int queryEnd,
  const char* pSubject, int subjectStart, int subjectEnd )
{

  outputStream_ << "Alignment score: " << alignment.totalScore_ << endl;

  int queryNext(queryStart), subjectNext(subjectStart);
  int* pProteinNext; 
  int* pDNANext;

  const char* pProteinChar;
  const char* pDNAChar;

  char queryChar, subjectChar, matchChar, codonChar;

  char* pProteinOut;
  char* pDNAOut;

  AlignBreakType proteinGap, DNAGap, DNAFrameShift;


  if (isQueryProtein_)
  {
    pProteinNext = &queryNext;
    pDNANext     = &subjectNext;
    pProteinChar = pQuery;
    pDNAChar     = pSubject;
    pProteinOut  = &queryChar;
    pDNAOut      = &subjectChar;
    proteinGap   = alignBreakGapQuery;
    DNAGap       = alignBreakGapSubject;
    DNAFrameShift = alignBreakFrameShiftSubject;
  } // ~if
  else
  {
    pProteinNext = &subjectNext;
    pDNANext     = &queryNext;
    pProteinChar = pSubject;
    pDNAChar     = pQuery;
    pProteinOut  = &subjectChar;
    pDNAOut      = &queryChar;
    proteinGap   = alignBreakGapSubject;
    DNAGap       = alignBreakGapQuery;
    DNAFrameShift = alignBreakFrameShiftQuery;
  } // ~else

  pCursorSeq1_ = &pBufSeq1_[numCols_]; // triggers reset to start of line

  *pBufSeq1_='\0'; // ensures 3 null strings are printed at first line break 

  Alignment::iterator i(alignment.begin());

  while( i!=alignment.end() ) 
  {
    //    cout << queryNext << " " << subjectNext << endl;
    
    if (i->numMatches>0)
    { 
      codonChar = getCodon( pDNAChar );
      matchChar = (codonChar==toupper(*pProteinChar)) ? '|' : 'x';
      codonChar = (matchChar=='|') ? '|' : codonChar;

      //      matchChar = (gResidueNames[ ttCodon[ ttDNA[ *(pDNAChar) ] << 4
      //                | ttDNA[ *(pDNAChar+1) ] << 2
      //                | ttDNA[ *(pDNAChar+2) ] ] ]
      //                   == toupper(*pProteinChar))
      //        ? '|' : 'x';


      *pProteinOut = *(pProteinChar++);
      *pDNAOut     = *(pDNAChar++);

      outputAlignmentColumn
        (queryChar, codonChar, subjectChar, queryNext, subjectNext ); 
      (*pProteinNext)++;
      (*pDNANext)++;

      *pProteinOut = '.';
      *pDNAOut     = *(pDNAChar++);
      outputAlignmentColumn
        (queryChar, matchChar, subjectChar, queryNext, subjectNext ); 
      (*pDNANext)++;

      *pDNAOut     = *(pDNAChar++);
      outputAlignmentColumn
        (queryChar, matchChar, subjectChar, queryNext, subjectNext ); 
      (*pDNANext)++;


      i->numMatches--;
    }
    else
    {
#ifdef XXX
      if (i->breakType==alignBreakMismatch)
      {
        //      outputAlignmentColumn
        //        (*(pQuery++), 'x', *(pSubject++), queryNext++, subjectNext++ ); 
      *pProteinOut = *(pProteinChar++);
      *pDNAOut     = *(pDNAChar++);
      outputAlignmentColumn
        (queryChar, 'x', subjectChar, queryNext++, subjectNext++ ); 

      *pProteinOut = '.';
      *pDNAOut     = *(pDNAChar++);
      outputAlignmentColumn
        (queryChar, 'x', subjectChar, queryNext++, subjectNext++ ); 

      *pDNAOut     = *(pDNAChar++);
      outputAlignmentColumn
        (queryChar, 'x', subjectChar, queryNext++, subjectNext++ ); 

      } // ~if
      else 
#endif
      if (i->breakType==DNAFrameShift)
      {
      *pProteinOut = '-';
      *pDNAOut     = *(pDNAChar++);
      outputAlignmentColumn
        (queryChar, ' ', subjectChar, queryNext, subjectNext ); 
      queryNext   += 1 * (!isQueryProtein_);
      subjectNext += 1 * ( isQueryProtein_);

      } // ~else if
      else if (i->breakType==proteinGap)
      {

        *pProteinOut = '-';
        *pDNAOut     = *(pDNAChar++);
        outputAlignmentColumn
          (queryChar, ' ', subjectChar, queryNext, subjectNext ); 
        queryNext   += 1 * (!isQueryProtein_);
        subjectNext += 1 * ( isQueryProtein_);

        *pDNAOut     = *(pDNAChar++);
        outputAlignmentColumn
          (queryChar, ' ', subjectChar, queryNext, subjectNext ); 
        queryNext   += 1 * (!isQueryProtein_);
        subjectNext += 1 * ( isQueryProtein_);

        *pDNAOut     = *(pDNAChar++);
        outputAlignmentColumn
          (queryChar, ' ', subjectChar, queryNext, subjectNext ); 
        queryNext   += 1 * (!isQueryProtein_);
        subjectNext += 1 * ( isQueryProtein_);

      } // ~else if
      else if (i->breakType==DNAGap)
      {

        *pProteinOut = *(pProteinChar++);
        *pDNAOut     = '-';
        
        outputAlignmentColumn
          (queryChar, ' ', subjectChar, queryNext, subjectNext ); 
        queryNext   += 1 * (isQueryProtein_);
        subjectNext += 1 * (!isQueryProtein_);

        *pProteinOut = '.';
        outputAlignmentColumn
          (queryChar, ' ', subjectChar, queryNext, subjectNext ); 
        outputAlignmentColumn
          (queryChar, ' ', subjectChar, queryNext, subjectNext ); 

      } // ~else if
      i++;
    } // ~else

  } // ~while

  *pCursorSeq1_='\0';
  *pCursorSeq2_='\0';
  *pCursorAlign_='\0';

  outputAlignmentLine();


} // ~void MatchAlignerTranslatedProtein::formatAlignment


void MatchAlignerTranslatedDNA::formatAlignment
( Alignment& alignment,
  const char* pQuery, int queryStart, int queryEnd,
  const char* pSubject, int subjectStart, int subjectEnd )
{

  outputStream_ << "Alignment score: " << alignment.totalScore_ << endl;

  int queryNext(queryStart);
  int subjectNext(subjectStart);
  char codonChar;

  pCursorSeq1_ = &pBufSeq1_[numCols_]; // triggers reset to start of line

  *pBufSeq1_='\0'; // ensures 3 null strings are printed at first line break 
  //  *pBufSeq2_='\0'; // - this line not necessary??
  //  *pBufAlign_='\0'; // - this line not necessary??

  Alignment::iterator i(alignment.begin());

  while( i!=alignment.end() ) 
  {
    //    cout << queryNext << " " << subjectNext << endl;
    
    if (i->numMatches>0)
    { // TBD do these en masse?
      codonChar = getCodon( pQuery );
      if ((codonChar!='X')&&(codonChar == getCodon(pSubject)))
      {
        outputAlignmentColumn
        (
         *(pQuery++),
         codonChar,
         *(pSubject++), 
         queryNext++, 
         subjectNext++ 
        ); 
        outputAlignmentColumn
          (*(pQuery++), '.', *(pSubject++), queryNext++, subjectNext++ ); 
        outputAlignmentColumn
        (*(pQuery++), '.', *(pSubject++), queryNext++, subjectNext++ ); 
      } // ~if
      else
      {
        outputAlignmentColumn
          ( *pQuery, matchChar(*pQuery,*pSubject) ? '|' : 'x', *pSubject, queryNext++, subjectNext++ );
        pQuery++; pSubject++;
        outputAlignmentColumn
          ( *pQuery, matchChar(*pQuery,*pSubject) ? '|' : 'x', *pSubject, queryNext++, subjectNext++ );
        pQuery++; pSubject++;
        outputAlignmentColumn
          ( *pQuery, matchChar(*pQuery,*pSubject) ? '|' : 'x', *pSubject, queryNext++, subjectNext++ );
        pQuery++; pSubject++;
      } // ~else

    i->numMatches--;
    }
    else
    {
#ifdef XXX
      if (i->breakType==alignBreakMismatch)
      {
        outputAlignmentColumn
          (*(pQuery++), 'x', *(pSubject++), queryNext++, subjectNext++ ); 
        outputAlignmentColumn
          (*(pQuery++), 'x', *(pSubject++), queryNext++, subjectNext++ ); 
        outputAlignmentColumn
          (*(pQuery++), 'x', *(pSubject++), queryNext++, subjectNext++ ); 
      } // ~if
      else 
#endif
      if (i->breakType==alignBreakFrameShiftQuery)
      {
        outputAlignmentColumn
          (*(pQuery++), ' ', '-', queryNext++, subjectNext ); 

      } // ~else if
      else if (i->breakType==alignBreakFrameShiftSubject)
      {
        outputAlignmentColumn
          ( '-', ' ', *(pSubject++), queryNext, subjectNext++ ); 
      } // ~else if
      i++;
    } // ~else

  } // ~while

  *pCursorSeq1_='\0';
  *pCursorSeq2_='\0';
  *pCursorAlign_='\0';

  outputAlignmentLine();



} // ~MatchAlignerTranslatedDNA::formatAlignment



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

void print( vector<ScoreType>& v )
{
  for( vector<ScoreType>::iterator i(v.begin());i!=v.end();++i)
    cout << "Score: " << *i << endl;

}


void print( PathMatrix<PathType>& p )
{
  for (int i(0); i< p.front().size(); i++)
  {
    for (vector<vector<PathType> >::iterator j(p.begin());
         j!=p.end();++j)
    {
      cout << (((*j)[i] & fromW)?'-':'.')
           << (((*j)[i] & fromSW)?'/':'.') 
           << (((*j)[i] & fromN)?'|':'.') 
           << "\t";
    } // ~for j   
    cout << endl;
  } // ~for i

} // ~print


ScoreType ColumnFillerBasic::operator()( PathMatrix<PathType>& matrix )
{
  ScoreType lastScore;
  int i,j;

  matrix.resize(bandLength_, vector<PathType>(colSize_) );

  // Fill in first column
  matrix[0][bandExtension_]=fromFinished;
  v1_[bandExtension_]=0;

  //  cout << "doing first column" << endl;
  for (j= bandExtension_+1; j< colSize_ ; j++)
  {
    v1_[j] = fillCell_( matrix[0][j],
                        veryBadScoreIndeed,
                        veryBadScoreIndeed,
                        v1_[j-1]+gapScore_);

  } // ~for i
  //  print (v1_);
  pLast_    = &v1_;
  pCurrent_ = &v2_;


  // Fill in next bandExtension_ columns

  for ( i=1 ; i <= bandExtension_ ; i++  ) // NB <= not <
  {
    //    cout << "doing leftmost columns" << i << endl;

    lastScore = fillCell_( matrix[i][bandExtension_-i],
                                 veryBadScoreIndeed,
                                 (*pLast_)[bandExtension_-i+1]+gapScore_,
                                 veryBadScoreIndeed );
    (*pCurrent_)[bandExtension_-i]=lastScore;
    for ( j=bandExtension_-i+1; j < colSize_-1 ; j++)
    {
      lastScore = fillCell_
        ( matrix[i][j],
          (*pLast_)[j]
          + getScore_(p1_[i-1],p2_[j-bandExtension_+i-1]),
          //      + ((tolower(p1_[i-1])==tolower(p2_[j-bandExtension_+i-1])) 
          //         ? matchScore : mismatchScore),
          (*pLast_)[j+1]    + gapScore_,
          lastScore + gapScore_ );
      (*pCurrent_)[j] = lastScore;
    } // ~for j

    (*pCurrent_)[j] = fillCell_
      ( matrix[i][j],
        (*pLast_)[j]
        + getScore_(p1_[i-1],p2_[j-bandExtension_+i-1]),
        // + ((tolower(p1_[i-1])==tolower(p2_[j-bandExtension_+i-1])) ?
        //   matchScore : mismatchScore),
        veryBadScoreIndeed,
        lastScore + gapScore_ );
    //  print (*pCurrent_);
    temp_ = pCurrent_; pCurrent_ = pLast_; pLast_ = temp_;
  } // ~for i


  // Fill in main set of columns
  for ( i=(bandExtension_+1) ; 
        i < bandLength_ -bandExtension_; 
        i++ ) 
  {
    //    cout << "doing centre columns" << i << endl;

    lastScore = veryBadScoreIndeed;

    for ( j=0; j < colSize_-1 ; j++)
    {
      lastScore = fillCell_
        ( matrix[i][j],
          (*pLast_)[j]
          + getScore_(p1_[i-1],p2_[j-bandExtension_+i-1]),
          //      + ((tolower(p1_[i-1])==tolower(p2_[j-bandExtension_+i-1])) 
          //   ? matchScore : mismatchScore),
          (*pLast_)[j+1]    + gapScore_,
          lastScore + gapScore_ );
      (*pCurrent_)[j]=lastScore;
    } // ~for j

    (*pCurrent_)[j] = fillCell_
      ( matrix[i][j],
        (*pLast_)[j]
        + getScore_(p1_[i-1],p2_[j-bandExtension_+i-1]),
        //      + ((tolower(p1_[i-1])==tolower(p2_[j-bandExtension_+i-1])) 
        // ? matchScore : mismatchScore),
          veryBadScoreIndeed,
          lastScore + gapScore_ );
    //  print (*pCurrent_);
    temp_ = pCurrent_; pCurrent_ = pLast_; pLast_ = temp_;

  } // ~for i

  // Fill in last bandExtension_ columns;


  for ( i=bandLength_-bandExtension_ ; i < bandLength_ ; i++ ) 
  {
    //    cout << "doing rightmost columns" << i << endl;

    lastScore = veryBadScoreIndeed;

    for ( j=0; j < colSize_-1-i-bandExtension_+bandLength_ ; j++)
    {
      lastScore = fillCell_
        ( matrix[i][j], 
          (*pLast_)[j]
          + getScore_(p1_[i-1],p2_[j-bandExtension_+i-1]),
          //      + ((tolower(p1_[i-1])
          //   ==tolower(p2_[j-bandExtension_+i-1])) 
          // ? matchScore : mismatchScore),
          (*pLast_)[j+1]    + gapScore_,
          lastScore + gapScore_ );
      (*pCurrent_)[j]=lastScore;
    } // ~for j
    //  print (*pCurrent_);
    temp_ = pCurrent_; pCurrent_ = pLast_; pLast_ = temp_;

  } // ~for i

  // Set position of the last cell filled in ready for
  // the traceback
  matrix.lastCell_.first
    = matrix.end()-1;
  matrix.lastCell_.second
    = matrix.lastCell_.first->begin() + (colSize_-bandExtension_-1);

  return (*pLast_)[colSize_-bandExtension_-1];

} // ~void ColumnFillerBasic::operator()


bool TraceBackerBasic::operator()
( PathMatrix<PathType>::CellIterator& current )
{
  //  cout << (int)(current.second-current.first->begin());

  alignment_.push_front( AlignInfo() );
  alignment_.front().breakType=alignBreakNull;

  if (*current.second&fromFinished)
  {
    //    cout << "finished traceback" << endl;
    return false;
  } // ~if
  else if (*current.second&fromW)
  {
    //    cout << "going W" << endl;
    current.second 
      = (current.first-1)->begin() 
      + (current.second-current.first->begin());
    --current.first;
    alignment_.front().numMatches++; // doesn't check for mismatch
  } // ~else if
  else if (*current.second&fromN)
  {
    //    cout << "going N" << endl;
    assert(current.second!=current.first->begin());
    --current.second;
    alignment_.push_front( AlignInfo() );
    alignment_.front().breakType=p1Gap_; 
  } // ~else if
  else if (*current.second&fromSW)
  {
    //    cout << "going SW" << endl;
    current.second 
      = (current.first-1)->begin() 
      + (current.second-current.first->begin());
    ++current.second;
    --current.first;
    alignment_.push_front( AlignInfo() );
    alignment_.front().breakType=p2Gap_; 
  } // ~else if
  else assert (1==0);

  return true;
} // ~bool TraceBackerBasic::operator()


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



  template<> ScoreType ColumnFiller3D::getScoreW<false>
  ( int j, int k, int l )
  {
    //    cout << "W blocked " << j << " " << k << " " << l << endl;
    return veryBadScoreIndeed;
  } // ~template<> ScoreType getScoreW<false>

  template<> ScoreType ColumnFiller3D::getScoreSW<false>
  ( int j, int k, int l )
  {
    //    cout << "SW blocked " << j << " " << k << " " << l << endl;
    return veryBadScoreIndeed;
  } // ~template<> ScoreType getScoreSW<false>

  template<> ScoreType ColumnFiller3D::getScoreN<false>
  ( int j, int k, int l )
  {
    //    cout << "N blocked " << j << " " << k << " " << l << endl;
    return veryBadScoreIndeed;
  } // ~template<> ScoreType getScoreN<false>

  template<> ScoreType ColumnFiller3D::getScoreChar<false>
  ( int i, int j, int k, int l)
  {
    //    cout << "char blocked " << i << " " << j << " " << k << " " 
    // << l << endl;
    return (ScoreType)0;
  }

// This specialisation fills in the first cell
  template<> void ColumnFiller3D::fillCell3D< false, false, false, false>
  ( PathMatrix<PathType3D>& matrix, int i, int j )
  {
    //    cout << "filling in first cell" << i << " " << j << endl;
    for (int k(0); k < numFrames1_; k++ )
    {
      for (int l(0); l < numFrames2_; l++ )
      {
        (*pCurrent_)[j][k][l] = (k+l)*frameShiftScoreBLOSUM;
        matrix[i][j][k][l] |= fromPrevFrame1 * (k!=0);
        matrix[i][j][k][l] |= fromPrevFrame2 * (l!=0);
        //              cout << i << " " << j << " " << k << " " << l << " "
        //         << (*pCurrent_)[j][k][l] << " "
        //          << (int)matrix[i][j][k][l] << endl;
      } // ~for l
    } // ~for k
    matrix[i][j][0][0] |= fromFinished;
  }

ColumnFiller3D::ColumnFiller3D
( const char** p1Trans, int p1Size, int p1FinalFrame,
  const char** p2Trans, int p2Size, int p2FinalFrame,
  int bandExtension,
  const ScoreTable& scoreTable ) :
  bandExtension_( (bandExtension<(p1Size/2))
                  ? bandExtension
                  : max( (p1Size/2)-1, 0 ) ), 
  //  bandExtension_(bandExtension), 
  bandWidth_(p2Size-p1Size+1),
  bandLength_(p1Size+1),
  colSize_(p2Size-p1Size+1+(2*bandExtension_)),
  fillCell_(),
  finalFrame1_(p1FinalFrame),
  finalFrame2_(p2FinalFrame),
  numFrames1_((p1Trans[1]==NULL)?1:gNumReadingFrames),
  numFrames2_((p2Trans[1]==NULL)?1:gNumReadingFrames),
  v1_(colSize_, veryBadScore3D ), 
  v2_(colSize_, veryBadScore3D ),
  pLast_(&v1_),
  pCurrent_(&v2_),
  getScore_(scoreTable)
{
  
  p1_[0] = p1Trans[0];
  p1_[1] = p1Trans[1];
  p1_[2] = p1Trans[2];
  p2_[0] = p2Trans[0];
  p2_[1] = p2Trans[1];
  p2_[2] = p2Trans[2];

  //    cout << "C3D::C3D " << bandExtension_ << " " << bandWidth_ << " "
  //   << bandLength_ << " " << colSize_ << " " 
  //    << numFrames1_ << " "
  //   << numFrames2_ << " "
  //   << finalFrame1_ << " "
  //   << finalFrame2_ << endl;

} // ~ColumnFiller3D::ColumnFiller3D







ScoreType ColumnFiller3D::operator()( PathMatrix<PathType3D>& matrix )
{

  //  ScoreType lastScore, prevFrameScore1, prevFrameScore2;
  int i,j,k,l;

  matrix.resize(bandLength_, vector<PathType3D>(colSize_) );

  //  matrix[0].resize(colSize_);
  
  pLast_    = &v2_;
  pCurrent_ = &v1_;

  // Fill in first column
  //  matrix[0][bandExtension_][0][0]=fromFinished;
  //  (*pCurrent_)[bandExtension_][0][0]=0;
  //  cout << "doing first cell" << endl;
  fillCell3D<false, false, false, false>( matrix, 0, bandExtension_);
  


  //    cout << "doing first column" << endl;
  for (j = (bandExtension_+1); j< colSize_ ; j++)
  {
    fillCell3D<false, false, true, false>( matrix, 0, j);
  } // ~for i

  temp_ = pCurrent_; pCurrent_ = pLast_; pLast_ = temp_;

  // Fill in next bandExtension_ columns

  for ( i=1 ; i <= bandExtension_ ; i++  ) // NB <= not <
  {
    //    cout << "doing leftmost columns" << i << endl;

    fillCell3D<false, true, false, false>( matrix, i, bandExtension_-i );


    for ( j=bandExtension_-i+1; j < colSize_-1 ; j++)
    {
      //      cout << "doing main chunk ";
      fillCell3D<true, true, true, true>( matrix, i, j);

    } // ~for j

    // Next if statement is unnecessary, as this scope is only executed
    // if band extension >= 1, since colSize > 2 * bandExtension
    //    if (colSize_>1)
    //   {
    //   cout << "doing last one ";
      fillCell3D<true, false, true, true>( matrix, i, j); // NB wot if loop never exec'd??
      //  }


    temp_ = pCurrent_; pCurrent_ = pLast_; pLast_ = temp_;
  } // ~for i


  // Fill in main set of columns
  if (colSize_==1)
  {
    for ( i=(bandExtension_+1) ; 
          i < bandLength_ -bandExtension_; 
          i++ ) 
    {
        fillCell3D<true, false, false, true>( matrix, i, 0);
        temp_ = pCurrent_; pCurrent_ = pLast_; pLast_ = temp_;
    } // ~for i
  } // ~if
  else
  {
    for ( i=(bandExtension_+1) ; 
          i < bandLength_ -bandExtension_; 
          i++ ) 
      {
        //      cout << "doing centre columns" << i << endl;

        //      cout << "doing first one ";
        fillCell3D<true, true, false, true>( matrix, i, 0);

        for ( j=1; j < colSize_-1 ; j++)
        {
          //      cout << "doing main chunk ";
          fillCell3D<true, true, true, true>( matrix, i, j);
        } // ~for j

        //      cout << "doing last one ";
        fillCell3D< true, false, true, true>(matrix, i, j);
        temp_ = pCurrent_; pCurrent_ = pLast_; pLast_ = temp_;
      } // ~for i
  } // ~else

  // Fill in last bandExtension_ columns;


  for ( i=bandLength_-bandExtension_ ; i < bandLength_ ; i++ ) 
  {
    //    cout << "doing rightmost columns" << i << endl;
    
    fillCell3D<true, true, false, true>( matrix, i, 0);

    for ( j=1; j < colSize_-1-i-bandExtension_+bandLength_ ; j++)
    {
      //      cout << "doing main chunk ";
      fillCell3D< true, true, true, true>( matrix, i, j);
    } // ~for j
    //    print (*pCurrent_);
    temp_ = pCurrent_; pCurrent_ = pLast_; pLast_ = temp_;

  } // ~for i

  // Set position of the last cell filled in ready for
  // the traceback
  matrix.lastCell_.first
    = matrix.end()-1;
  matrix.lastCell_.second
    = matrix.lastCell_.first->begin()+ (colSize_-bandExtension_-1);

  return (*pLast_)[colSize_-bandExtension_-1][finalFrame1_][finalFrame2_];

} // ~void ColumnFiller3D::operator()


bool TraceBacker3D::operator()
( PathMatrix<PathType3D>::CellIterator& current )
{

  //    cout << (int)(current.second-current.first->begin())
  // << " frame1: " << k_ << " frame2: " << l_ 
  // << endl;

  alignment_.push_front( AlignInfo() );
  alignment_.front().breakType=alignBreakNull;

  if (((*current.second)[k_][l_])&fromFinished)
  {
    return false;
  }
  else if (((*current.second)[k_][l_])&fromW)
  {
    //  cout << "going W" << endl;
    current.second 
      = (current.first-1)->begin() 
      + (current.second-current.first->begin());
    --current.first;
    alignment_.front().numMatches++;
  } // ~else if
  else if (((*current.second)[k_][l_])&fromN)
  {
    //    cout << "going N" << endl;
    assert(current.second!=current.first->begin());
    --current.second;
    alignment_.push_front( AlignInfo() );
    alignment_.front().breakType=p1Gap_;
  } // ~else if
  else if (((*current.second)[k_][l_])&fromSW)
  {
    //    cout << "going SW" << endl;
    current.second 
      = (current.first-1)->begin() 
      + (current.second-current.first->begin());
    ++current.second;
    --current.first;
    alignment_.push_front( AlignInfo() );
    alignment_.front().breakType=p2Gap_;
  } // ~else if
  else if (((*current.second)[k_][l_])&fromPrevFrame1)
  {
    //    cout << "shifting frame p1" << endl;
    if (k_==0)
    { // then need to do a SW shift
      // (equivalent to a W shift in the full d.p. matrix)
      current.second 
        = (current.first-1)->begin() 
        + (current.second-current.first->begin());
      ++current.second;
      --current.first;
      k_=2;
    } // ~if
    else k_--;
    alignment_.push_front( AlignInfo() );
    alignment_.front().breakType=p1FrameShift_;
  } // ~else if
  else if (((*current.second)[k_][l_])&fromPrevFrame2)
  {
    //    cout << "shifting frame p2" << endl;
    if (l_==0)
    { // then need to do a N shift
      // (equivalent to a N shift in the full d.p. matrix)
      assert(current.second!=current.first->begin());
      --current.second;
      l_=2;
    } // ~if
    else l_--;
    alignment_.push_front( AlignInfo() );
    alignment_.front().breakType=p2FrameShift_;
  } // ~else if
  else assert (1==0);

  return true;

} // ~bool TraceBacker3D::operator()




// End of file MatchAligner.cpp

---