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

// Description:

// Includes:

#include "MatchStore.h"
#include "HashTableGeneric.h"
#include "QueryManager.h"
#include "MatchStoreGapped.h"
#include "SequenceReader.h"
#include <iomanip> // for setprecision() etc...

// Match member function definitions


const char* MatchImp::getSubjectName( void ) const
{
  //  cout << "names_" << myStore_->names_.size() << endl; // %%%
  //  assert(myStore_->names_.find(subjectNum_)!=myStore_->names_.end());
  //  return myStore_->names_[ subjectNum_ ];

  return myStore_->reader_.getSequenceName(subjectNum_);

  //  string name;
  //  myStore_->subjectTable_.getSequenceName( name, subjectNum_ );
  //  return name;
}

SequenceNumber MatchImp::getQueryNum( void ) const
{
  return myStore_->queryNum_;
}

string MatchImp::getQueryName( void ) const
{
  return myStore_->queryName_;
}

int  MatchImp::getQuerySize( void ) const
{
  return myStore_->queryBases_;
}

// MatchStore member functions definitions


  // Function Name: printResult
  // Arguments: const QueryResult& (in), ostream& (out)
  // Convert the contents of result to ASCII and send to the specified
  // output stream
void MatchStoreImp::printResult
( ostream& outputStream  ) const
{ 
  
  if (empty()) return;

  outputStream << endl << "Matches For Query " 
               << queryNum_
               << " (" << queryBases_
               << " bases): " << queryName_
               << "\n\n";    
 
  for( vector<MatchImp*>::const_iterator i(imps_.begin());i!=imps_.end();++i)
  {
    //    MatchImp* i(*j);

    outputStream << (((*i)->isQueryForward_)?"F":"R")
                 <<" " << (*i)->subjectNum_
                 << "\t: " << (*i)->getSubjectName()// names_[(*i)->subjectNum_]
                 << "\tBases: " << (*i)->numBases_
                 << "\tQ: " << (*i)->queryStart_ 
                 << " to " << (*i)->queryEnd_
                 << "\tS: " << (*i)->subjectStart_ 
                 << " to " << (*i)->subjectEnd_ 
                 << "\n"; 

  }  // ~for i
   
  outputStream << endl;

} // ~QueryManager::printResult

MatchStoreImp::~MatchStoreImp()
{
  for (vector<MatchImp*>::iterator i(imps_.begin()) ; i!= imps_.end() ; ++i )
    delete (*i);

}

void MatchStoreImp::clear( void )
{ 
    
  for (vector<MatchImp*>::iterator i(imps_.begin()) ; i!= imps_.end() ; ++i )
    delete (*i);
  vector<Match*>::clear();
  imps_.clear();
  //    matchImps_.clear();
  queryBases_ = 0;
}




void MatchStoreImp::addMatch( 
                          SequenceNumber subjectNum,
                          SequenceOffset numBases,
                          SequenceOffset queryStart,
                          SequenceOffset queryEnd,
                          SequenceOffset subjectStart,
                          SequenceOffset subjectEnd,
                          bool isQueryForward,
                          bool isSubjectForward )
{
  DEBUG_L2(             subjectNum << " " << 
                        numBases << " " << 
                        queryStart << " " << 
                        queryEnd << " " << 
                        subjectStart << " " << 
                        subjectEnd << " " << 
                        isQueryForward );

  imps_.push_back( new MatchImp ( this,
                                  subjectNum, 
                                  numBases, 
                                  queryStart, 
                                  queryEnd, 
                                  subjectStart, 
                                  subjectEnd, 
                                  isQueryForward,
                                  isSubjectForward) );
  //  push_back(static_cast<Match*>(imps_.back()));
  //  names_.insert(pair<SequenceNumber, string>(subjectNum, subjectName));
  //  cout << "inserting subject num" << subjectNum << " " << subjectName
  //     << endl;
  //  push_back(&matchImps_.back());
  //  MatchStore::iterator i(&back());
  //  DEBUG_L1( matchImps_.size() << " " << size() << " " 
  //    << i << " " << *i << " "  );
  //  DEBUG_L1( (*i)->getQueryNum() << " " << (*i)->getQueryName() << " " 
  //    << (*i)->getQueryStart() << " " << (*i)->getQueryEnd() << " " 
  //    << (*i)->getSubjectNum() << " " << (*i)->getSubjectName() << " " 
  //    << (*i)->getSubjectStart() << " " << (*i)->getSubjectEnd() ); 
  

}

void MatchStoreImp::setup( void )
{
  vector<Match*>::clear();
  for ( vector<MatchImp*>::iterator i(imps_.begin());i!=imps_.end();++i)
    push_back( static_cast<Match*>(*i));
}


// MatchTask member function definitions


void MatchTaskPrint::operator()( MatchStore& store )
{
  if (store.empty()) return;

  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);

  for( ;i!=store.end();++i)
  {

    outputStream_ << (((*i)->isQueryForward() )?"F":"R")
                  << (((*i)->isSubjectForward() )?"F":"R")
                  << " " << (*i)->getSubjectNum() 
                  << "\t: " << (*i)->getSubjectName()
                  << "\tScore: " << (*i)->getNumBases()
                  << "\tQ: " << (*i)->getQueryStart()  
                  << " to " << (*i)->getQueryEnd() 
                  << "\tS: " << (*i)->getSubjectStart()  
                  << " to " << (*i)->getSubjectEnd()  
                  << "\t" << 100.0*((*i)->getNumBases()) / 
      ((*i)->getQueryEnd()-(*i)->getQueryStart()+1)
                  << "\%\n"; 

  }  // ~for i
   
  outputStream_ << endl;


}

void MatchTaskPrintTabbed::operator()( MatchStore& store )
{
  if (store.empty()) return;

  outputStream_ << setprecision(2) << setiosflags(ios::fixed);
  for (MatchStore::iterator i(store.begin()); i!=store.end() ; ++i )
  outputStream_ 
    << ((*i)->isQueryForward()?"F":"R")
    << ((*i)->isSubjectForward()?"F":"R") << "\t"
    << (*i)->getQueryName() << "\t"
    << (*i)->getQueryStart() << "\t" 
    << (*i)->getQueryEnd() << "\t"
    << (*i)->getSubjectName() << "\t"
    << (*i)->getSubjectStart() << "\t"
    << (*i)->getSubjectEnd() << "\t"
    << (*i)->getNumBases() << "\t" 
    << 100.0*((*i)->getNumBases()) / 
       ((*i)->getQueryEnd()-(*i)->getQueryStart()+1)
    << endl;

}
void MatchTaskPrintReverse::operator()( MatchStore& store )
{
  if (store.empty()) return;

  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);

  for( ;i!=store.end();++i)
  {
    outputStream_ << (((*i)->isQueryForward() )?"F":"R")
                  << (((*i)->isSubjectForward() )?"F":"R")
                  <<" " << (*i)->getSubjectNum() 
                  << "\t: " << (*i)->getSubjectName()
                  << "\tScore: " << (*i)->getNumBases()
                  << "\tQ: " 
                  << ( (*i)->isQueryForward() ? (*i)->getQueryStart()  
                     : (*i)->getQuerySize()-(*i)->getQueryEnd()+1 )     
                  << " to " 
                  << ( (*i)->isQueryForward() ? (*i)->getQueryEnd()  
                     : (*i)->getQuerySize()-(*i)->getQueryStart()+1 )     
                  << "\tS: " << (*i)->getSubjectStart()  
                  << " to " << (*i)->getSubjectEnd()  
                  << "\t" << 100.0*((*i)->getNumBases()) / 
                    ((*i)->getQueryEnd()-(*i)->getQueryStart()+1)
                  << "\%\n"; 

  }  // ~for i
   
  outputStream_ << endl;


}

void MatchTaskPrintTabbedReverse::operator()( MatchStore& store )
{
  if (store.empty()) return;
  outputStream_ << setprecision(2) << setiosflags(ios::fixed);
  for (MatchStore::iterator i(store.begin()); i!=store.end() ; ++i )
  outputStream_ 
    << ((*i)->isQueryForward()?"F":"R") 
    << ((*i)->isSubjectForward()?"F":"R") << "\t"
    << (*i)->getQueryName() << "\t"
    << (   (*i)->isQueryForward() ? (*i)->getQueryStart()  
         : (*i)->getQuerySize()-(*i)->getQueryEnd()+1 ) << "\t"     
    << (   (*i)->isQueryForward() ? (*i)->getQueryEnd()  
         : (*i)->getQuerySize()-(*i)->getQueryStart()+1 ) << "\t"     
    << (*i)->getSubjectName() << "\t"
    << (*i)->getSubjectStart() << "\t"
    << (*i)->getSubjectEnd() << "\t"
    << (*i)->getNumBases() << "\t" 
    << 100.0*((*i)->getNumBases()) / 
       ((*i)->getQueryEnd()-(*i)->getQueryStart()+1)
    << endl;

}

#ifdef MOVED_TO_ALIGNERDOTCPP
// Definitions for MatchTaskAlign


PathType fromW = (PathType)1;
PathType fromN = (PathType)2;
PathType fromSW = (PathType)4;
PathType fromMismatch = (PathType)8;

// Next two definitions mean that the absolute value of each
// score in a score table need to be either less than 1024 
// or a power or two
ScoreType veryBadScoreIndeed = -67108864;
//ScoreType matchMask = 1024; // set this bit if characters match

ScoreType matchScore = 1; // | matchMask;
ScoreType mismatchScore= -2;
ScoreType gapStartScore = -4;
ScoreType gapContScore = -3;
ScoreType nScore = 0; // same as crossmatch defaults
//ScoreType veryBadScoreIndeed = -100000000;




const AlignBreakType alignBreakNull = 0;
const AlignBreakType alignBreakMismatch = 1;
const AlignBreakType alignBreakGapQuery = 2;
const AlignBreakType alignBreakGapSubject = 3;

MatchTaskAlign::MatchTaskAlign
( SourceReader& querySource, 
  SourceReader& subjectSource,
  int numCols,
  ostream& outputStream, 
  bool reverseQueryCoords,
  bool doAlignment ):
  querySource_(querySource), 
  subjectSource_(subjectSource),
  numCols_(numCols),
  reverseQueryCoords_(reverseQueryCoords),
  doAlignment_(doAlignment),
  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';
} // ~MatchTaskAlign::MatchTaskAlign


MatchTaskAlign::~MatchTaskAlign()
{
  delete [] pBufSeq1;
  delete [] pBufSeq2;
  delete [] pBufAlign;
} // ~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;
    subjectSource_.extractSource( &pSubject, //subjectData, 
                                  (*i)->getSubjectNum(),
                                  (*i)->getSubjectStart(),
                                  (*i)->getSubjectEnd() );

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

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

    align( 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)


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

  deque<AlignInfo> alignment;
  createAlignment
    ( alignment, 
      pQuery, queryStart, queryEnd, 
      pSubject, subjectStart, subjectEnd );
  formatAlignment
    ( alignment, 
      pQuery, queryStart, queryEnd, 
      pSubject, subjectStart, subjectEnd );

} // ~void MatchTaskAlign::align




void MatchTaskAlign::createAlignment
( 
 deque<AlignInfo>& 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;

  AlignBreakType gapInLargest, gapInSmallest;

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

  //  const char* p = pQuery;
  //  for (int i(queryStart); i<= queryEnd ; i++ ) cout << *p++;
  //  cout << endl;
  //  p = pSubject;
  //  for (int i(subjectStart); i<= subjectEnd ; i++ ) cout << *p++;
  // cout << endl;

  //  fillScoreTable(table);  // TBD do this in constructor

  vector<vector<PathType> > path;

  // set up data structure to store path
  path.resize( bandWidth+1 );

  // zeroth vector is kept null, to make indexing easier
  for (int i(1) ; i<=bandWidth; i++ ) path[i].resize( bandLength );

  // set up storage for current and previous column of scores
  vector<ScoreType> scoresCurrent, scoresLast;
  scoresCurrent.resize(bandWidth+2); // %%
  scoresLast.resize(bandWidth+2); // %%%

  // These stop the traceback going 'off track'
  scoresCurrent.front()=veryBadScoreIndeed;
  scoresCurrent.back()=veryBadScoreIndeed;

  scoresLast.front()=veryBadScoreIndeed;
  scoresLast.back()=veryBadScoreIndeed;


  vector<ScoreType>* pLast(&scoresLast);
  vector<ScoreType>* pCurrent(&scoresCurrent);
  vector<ScoreType>* pTemp;

  // set up storage to decide maximum score

  vector<pair<ScoreType, PathType> > scoresPossible;

  //  scoresPossible.resize(4); 
  scoresPossible.resize(3); 

  int i,j;

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

  ScoreType bestScore, thisScore;
  bool isMatch;

  for (i=0; i<bandLength; i++)
  {
    for (j=1 ; j <= bandWidth ; j++ ) // %%
    {
      //            cout << "j: " << j << " i: " << i << endl;
      path[j][i]=0; // TBD probably not needed

      // fill in direction indicators
      //    scoresPossible[0].second=0;

      scoresPossible[0].second=fromW;
      scoresPossible[1].second=fromN;
      scoresPossible[2].second=fromSW;

      // fill in corresponding scores

      //      scoresPossible[0].first=0;

      //      cout << p2[i+j-1] << p1[i] << endl;

      a[0]=p2[i+j-1];
      a[1]='-';
      scoresPossible[1].first=(*pCurrent)[j-1]+(*pTable_)[*b];

      a[1]=p1[i];
      //      isMatch=(*pTable_)[*b]&matchMask;
      isMatch=(tolower(a[0])==tolower(a[1]));

      //      scoresPossible[0].first=(*pLast)[j]+(*pTable_)[*b]-(isMatch*matchMask);
      scoresPossible[0].first=(*pLast)[j]+(*pTable_)[*b];
      //  cout << a[0] << a[1] << " " << (*pTable_)[*b]  << endl;
      // flag mismatches
      //      path[j][i] |= ((*pTable_)[*b]!=matchScore)*fromMismatch;
      path[j][i] |= (!isMatch)*fromMismatch;

      a[0]='-';

      scoresPossible[2].first=(*pLast)[j+1]+(*pTable_)[*b];

      //        cout << "scores 0/W/N/SW: " << scoresPossible[0].first << " " 
      //    << scoresPossible[1].first << " " 
      //    << scoresPossible[2].first << endl;

      sort(scoresPossible.begin(),scoresPossible.end());

      bestScore = scoresPossible.back().first;

      //      cout << "best score: " << bestScore << endl;

      (*pCurrent)[j]=bestScore;


      path[j][i] |=     
        scoresPossible[2].second; // * (bestScore!=0);

      path[j][i] |=  
        scoresPossible[1].second 
        * (scoresPossible[1].first==bestScore);
        //      * (bestScore!=0);

      path[j][i] |= 
        scoresPossible[0].second 
        * (scoresPossible[0].first==bestScore);
        //      * (bestScore!=0);

      //      path[j][i] |= 
      //        scoresPossible[0].second 
      //        * (scoresPossible[0].first==bestScore)
      //        * (bestScore!=0);

      //      cout << "path entry: " << (int)path[j][i] << endl;

    } // ~for j
    pTemp=pLast;
    pLast=pCurrent;
    pCurrent=pTemp;

  } // ~for i

  //    cout << bandWidth << " " << bandLength << endl;

     for (j=1 ; j <= bandWidth ; j++ )
     {
       //    cout << "j: ";
     for (i=0; i<bandLength; i++)
     {
       //       cout << (int)path[j][i] << " - ";
     } // ~for i
     //   cout << endl;
    } // ~for j


     cout << "alignment score: " << (*pLast)[bandWidth] << endl; 

  // et finalement, le traceback

  i= bandLength-1; j=bandWidth;

  alignment.push_front();
  alignment.front().breakType=alignBreakNull;

  while ((i>=0)&&(j>=1))
  {
    //      cout << j << " " << i << "... " << endl;
    if ((path[j][i]&fromW)&&(i>=0))
    {
      if (path[j][i]&fromMismatch)
      {
        //      cout << " mismatch" << endl;
        alignment.push_front();
        alignment.front().breakType=alignBreakMismatch;
      } // ~if
      else 
      {
        //      cout << " match" << endl;
        alignment.front().numMatches++;
      } // ~else
      i--;
    } // ~if
    else if ((path[j][i]&fromN)&&(j>=1))
    {
      // cout << " gap in smallest" << endl;
      j--;
      alignment.push_front();
      alignment.front().breakType=gapInSmallest; 
    } // ~else if
    else if ((path[j][i]&fromSW)&&(j<bandWidth)&&(i>0))
    {
      //    cout << " gap in largest" << endl;
      i--; j++;
      alignment.push_front();
      alignment.front().breakType=gapInLargest;
    } // ~else if
    //    else assert(1==0) //if (i>0)
      //    {
      //   cout << " mismatch" << endl;
      //  i--;
      //  alignment.push_front();
      // alignment.front().breakType=alignBreakMismatch;
      // }
    else assert(1==0); // stuck! - shouldn't happen!

  } // ~while
  //  cout << "final i,j=" << i << ", " << j << endl;
  for ( ; j > 1 ; j-- )
  {
      alignment.push_front();
      alignment.front().breakType=gapInSmallest; 
  }
} // ~void MatchTaskAlign::createAlignment


void MatchTaskAlign::formatAlignment
( 
 deque<AlignInfo>& alignment,
 const char* pQuery, int queryStart, int queryEnd,
 const char* pSubject, int subjectStart, int subjectEnd )
{
  //  cout << "format alignment" << endl;

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

  char* pCursorSeq1(&pBufSeq1[numCols_]); // triggers reset to start of line
  char* pCursorAlign(NULL);
  char* pCursorSeq2(NULL);

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

  deque<AlignInfo>::iterator i(alignment.begin());

  while( i!=alignment.end() ) 
  {
    //    cout << queryNext << " " << subjectNext << endl;
    if (pCursorSeq1==&pBufSeq1[numCols_])
    {
      //      cout << "about to write line" << endl;
      if (*pBufSeq1!='\0')
        outputAlignmentLine();
        //      cout << pBufSeq1 << endl << pBufAlign << endl << pBufSeq2 << endl<< endl;
      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
    
    if (i->numMatches>0)
    { // TBD do these en masse using memcpy 
    *(pCursorSeq1++)=*(pQuery++);
    *(pCursorAlign++)='|';
    *(pCursorSeq2++)=*(pSubject++);
    queryNext++;
    subjectNext++;
    i->numMatches--;
    }
    else
    {
      if (i->breakType==alignBreakMismatch)
      {
        *(pCursorSeq1++)=*(pQuery++);
        *(pCursorAlign++)='x';
        *(pCursorSeq2++)=*(pSubject++);
        queryNext++;
        subjectNext++;
      } // ~if
      else if (i->breakType==alignBreakGapQuery)
      {
        *(pCursorSeq1++)='-';
        *(pCursorAlign++)=' ';
        *(pCursorSeq2++)=*(pSubject++);
        //      queryNext++;
        subjectNext++;
      } // ~else if
      else if (i->breakType==alignBreakGapSubject)
      {
        *(pCursorSeq1++)=*(pQuery++);
        *(pCursorAlign++)=' ';
        *(pCursorSeq2++)='-';
        queryNext++;
        //      subjectNext++;
      } // ~else if
      i++;
    } // ~else

  } // ~while

  // check for unfinished lines
  //  if (pCursorSeq1!=&pBufSeq1[numCols_-1])
  //  {
  //  *pCursorSeq1++='\n';
  //  *pCursorSeq2++='\n';
  //  *pCursorAlign++='\n';
  *pCursorSeq1='\0';
  *pCursorSeq2='\0';
  *pCursorAlign='\0';
  //  }
  //  if (pCursorSeq1!=&pBufSeq1[numCols_])
  //  {
  //  }
  //  cout << pBufSeq1 << endl << pBufAlign << endl << pBufSeq2 << endl << endl;
  outputAlignmentLine();

} // ~void MatchTaskAlign::formatAlignment

void MatchTaskAlign::outputAlignmentLine( void )
{
  outputStream_ << pBufSeq1 << endl << pBufAlign << endl 
                << pBufSeq2 << endl << endl;
} // ~void MatchTaskAlign::outputAlignmentLine( void )



// Definitions for MatchTaskAlignDNA

bool MatchTaskAlignDNA::isTableFilled_ = false;
ScoreTable MatchTaskAlignDNA::table_;


MatchTaskAlignDNA::MatchTaskAlignDNA
( SourceReader& querySource, 
  SourceReader& subjectSource,
  int numCols=80,
  ostream& outputStream=cout,
  bool reverseQueryCoords=true,
  bool doAlignment=true ) : 
  MatchTaskAlign( querySource, 
                  subjectSource, 
                  numCols, 
                  outputStream,
                  reverseQueryCoords,
                  doAlignment )
{
  pTable_ = &table_;
  if (!isTableFilled_)
  {
    fillScoreTable();
    isTableFilled_=true;
  }

} // ~MatchTaskAlignDNA::MatchTaskAlignDNA


void MatchTaskAlignDNA::fillScoreTable( void )
{

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

  // TBD PAM/BLOSUM for proteins

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

  ScoreType thisScore;
  for (int i(0); i < 4; i++ )
  {
    for (int j(0); j < 4; j++ )
    {
      thisScore = ((i==j) ? matchScore: mismatchScore);
      a[0]=base[i]; a[1] = base[j]; table_[ *b ] = thisScore;
      //      a[0]=base[i]; a[1] = toupper(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[0]=toupper(base[i]); a[1] = toupper(base[j]); table_[ *b ] = thisScore;
      a[1] = toupper(base[j]); table_[ *b ] = thisScore;
    } // ~for j

    a[0]=base[i];  
    a[1] = 'n'; table_[ *b ] = nScore;
    a[1] = 'N'; table_[ *b ] = nScore;
    a[1] = '-'; table_[ *b ] = gapStartScore;

    a[0]=toupper(base[i]);  
    a[1] = 'n'; table_[ *b ] = nScore;
    a[1] = 'N'; table_[ *b ] = nScore;
    a[1] = '-'; table_[ *b ] = gapStartScore;

    a[1]=base[i];  
    a[0] = 'n'; table_[ *b ] = nScore;
    a[0] = 'N'; table_[ *b ] = nScore;
    a[0] = '-'; table_[ *b ] = gapStartScore;

    a[1]=toupper(base[i]);  
    a[0] = 'n'; table_[ *b ] = nScore;
    a[0] = 'N'; table_[ *b ] = nScore;
    a[0] = '-'; table_[ *b ] = gapStartScore;
    


  } // ~for i;

    a[0]='n'; a[1] = 'n'; table_[ *b ] = nScore;
    a[0]='n'; a[1] = 'N'; table_[ *b ] = nScore;
    a[0]='N'; a[1] = 'n'; table_[ *b ] = nScore;
    a[0]='N'; a[1] = 'N'; table_[ *b ] = nScore;

    a[0]='-'; a[1] = 'n'; table_[ *b ] = gapStartScore;
    a[0]='-'; a[1] = 'N'; table_[ *b ] = gapStartScore;
    a[0]='n'; a[1] = '-'; table_[ *b ] = gapStartScore;
    a[0]='N'; a[1] = '-'; table_[ *b ] = gapStartScore;

    a[0]='-'; a[1] = '-'; table_[ *b ] = gapContScore;

} // ~void MatchTaskAlignDNA::fillScoreTable( void )

// Definitions for MatchTaskAlignProtein
 
bool MatchTaskAlignProtein::isTableFilled_ = false;
ScoreTable MatchTaskAlignProtein::table_;

MatchTaskAlignProtein::MatchTaskAlignProtein
( SourceReader& querySource, 
  SourceReader& subjectSource,
  int numCols=80,
  ostream& outputStream=cout ) : 
  MatchTaskAlign( querySource, subjectSource, numCols, outputStream, false, true )
{
  pTable_ = &table_;
  if (!isTableFilled_)
  {
    fillScoreTable();
    isTableFilled_=true;
  }

} // ~MatchTaskAlignProtein::MatchTaskAlignProtein

ScoreType blosumScoreStopCodon= -4; // as used by NCBI - see URL in .h file
ScoreType blosumScoreGapStart= -10;
ScoreType blosumScoreGapExtend=-1;


void MatchTaskAlignProtein::fillScoreTable( void )
{

  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 ] = blosumScoreStopCodon; // a0=uc, a1=*
    a[1] = '-'; table_[ *b ] = blosumScoreGapStart;  // 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 ] = blosumScoreStopCodon; // a0=lc, a1=*
    a[1] = '-'; table_[ *b ] = blosumScoreGapStart;  // a0=lc, a1=-

    a[1] = blosumResidues[i];
    a[0] = '*'; table_[ *b ] = blosumScoreStopCodon; // a0=*, a1=uc
    a[0] = '-'; table_[ *b ] = blosumScoreGapStart;  // a0=-, a1=uc

    a[1] = tolower(a[1]);
    a[0] = '*'; table_[ *b ] = blosumScoreStopCodon; // a0=*, a1=lc
    a[0] = '-'; table_[ *b ] = blosumScoreGapStart;  // a0=-, a1=lc

    ++pScore;
  } // ~for i

  a[0] = '*'; a[1] = '*'; table_[ *b ] = 1;
  a[0] = '*'; a[1] = '-'; table_[ *b ] = blosumScoreGapStart;
  a[0] = '-'; a[1] = '*'; table_[ *b ] = blosumScoreGapStart;
  a[0] = '-'; a[1] = '-'; table_[ *b ] = blosumScoreGapExtend;

} // ~void MatchTaskAlignProtein::fillScoreTable( void )
#endif

// ### Function Definitions ###

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

// End of file MatchStore.cpp

---