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

// Description:

// Includes:

#include "MatchStoreGapped.h"
#include "SequenceReader.h"
#include "HashTableGeneric.h"
#include "HashTablePacked.h"
#include "GlobalDefinitions.h"
#include <algorithm>

// ### Function Definitions ###

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

// MatchAlgorithm member function definitions

void MatchAlgorithm::operator()  
  ( WordSequence& querySeq, 
    MatchAdder& addMatch, 
    HashTableGeneric& subjectTable,
    int wordLength,
    int stepLength )
{
  //  cout << "MatchAlgorithm()" << endl;
  HitListVector hitList;
  wordLength_ = ((wordLength==-1)?subjectTable.getWordLength():wordLength);
  stepLength_ = ((stepLength==-1)?wordLength_:stepLength);
  assert(hitList.size()==0);

  if (dynamic_cast<HashTablePacked*>(&subjectTable)!=NULL)
  {
    sortNeeded_ = false;
  } // ~if

  generateHits( querySeq, hitList, subjectTable );
  //  cout << "MatchAlgorithm() - hits got: " << hitList.size() << endl;
  generateMatches( hitList, addMatch );
}

void MatchAlgorithm::generateHits
( WordSequence& querySeq, HitListVector& hitList, HashTableGeneric& subjectTable )
{
  subjectTable.setNumRepeats( numRepeats_ );
  //  wordLength_ = subjectTable.getWordLength();
  //  stepLength_ = subjectTable.getStepLength();
  subjectTable.matchSequence( querySeq, hitList );
}



void MatchAlgorithmGapped::generateMatches
  ( HitListVector& hitList, MatchAdder& addMatch ) 
{

  if (hitList.size()==0) return;
  
  if (sortNeeded_)
  {
    sort( hitList.begin(), hitList.end(), LessThanSubject() );
  } // ~if 

  minHitsForMatch_ = max(2,minToProcess_ / stepLength_);
  maxQueryDiff_ = maxGap_ + stepLength_;

  bool inRun(false);

  HitListVector::iterator pFirstMatch( hitList.begin() );
  
  for ( int i(0) ; i < ((int)hitList.size()) - 1 ; ++i )
  {
    if (    ( hitList[i+1].subjectNum == hitList[i].subjectNum ) 
              && ( hitList[i+1].diff - hitList[i].diff  <=maxInsert_ )) 
    {
      if ( inRun==false ) 
      { 
        pFirstMatch = static_cast<HitList::iterator>(&hitList[i]); 
        inRun=true;
      } // ~if
    } // ~if
    else 
    {
      if (inRun==true)
        findMatchesInRange
        ( pFirstMatch, 
          static_cast<HitList::iterator>(&hitList[i+1]),
          addMatch );
      inRun=false;
    } // ~else if

  } // ~for

  if (inRun==true)
    findMatchesInRange( pFirstMatch, hitList.end(), addMatch );

  return;

} // ~findMatch


void MatchAlgorithmGapped::findMatchesInRange
( HitListVector::iterator first, HitListVector::iterator last, 
  MatchAdder& addMatch )
{

  if (first == last) return;
  // This means that single hit matches will not be reported, but you
  // can get these with MatchStoreUngapped

  sort( first, last, LessThanQuery() );

  //  cout << "sortedHits:" <<endl;
  //  for (HitListVector::iterator j(first); j!=last ; ++j)
  // cout << j << " " << j->subjectNum << ": " << j->queryPos << " " 
  // << j->diff+j->queryPos << endl;


  int lastQueryPos, numBases, gapSize; 
  SequenceOffset subjectStart, subjectEnd;

  while ( last-first >= minHitsForMatch_ )
  {
    //    cout << last-first << " left to go.\n";
    lastQueryPos = first->queryPos;
    numBases = wordLength_;
    HitListVector::iterator i(first);

    while(++i!=last)
    {

      gapSize = i->queryPos - lastQueryPos;
      //   cout << lastQueryPos << "|" << i->queryPos << endl; 
      if ( gapSize > maxQueryDiff_ ) break;
      numBases += min (wordLength_, gapSize );
      lastQueryPos = i->queryPos;
    } // ~while
  
    if (numBases>=minToProcess_)
    { 
      i--;
      //      cout << "addMatch: "  <<     first->queryPos << " " <<
      //  i->queryPos + wordLength_ - 1 << " " <<
      //  first->diff + first->queryPos << " " <<
      //i->diff + i->queryPos + wordLength_ - 1 << endl;

      subjectStart = first->diff + first->queryPos;
      subjectEnd= i->diff + i->queryPos + wordLength_ - 1; 

      // Put in the test below, as repetitive queries can occasionally
      // cause it to be false, which causes a crash at the alignment stage
      // TC 29.5.2
      if (subjectEnd>subjectStart)
      {
        addMatch 
        ( 
         first->subjectNum,
         numBases,
         first->queryPos,
         i->queryPos + wordLength_ - 1,
         subjectStart,
         subjectEnd
         );
      } // ~if
      i++;
    } // ~if

    first = i;
  
  } // ~while

} // ~findMatchesInRange





  // Function Name: findMatch
  // Arguments: vector<Match>& (out), int (in), int (in), int (in)
  // Returns: true if more than minToPrint bases were found in total.
  // Sorts its member HitInfo instances in the following way:
  // 1 (highest priority) by subjectNum: thus the hits for each subject
  // sequence are grouped together
  // 2 by diff: for a given subject sequence, hits which have the same 'diff' 
  // value are part of the same ungapped alignment between the subject
  // sequence and the query sequence.
  // 3 (lowest priority) by queryPos
  // After sorting, the routine looks for a region of hits which have the same
  // subjectNum and diff values but for which the queryPos value differs
  // by stepLength. Each such region denotes a consecutive series of matching 
  // base pairs between the query and subject, and the details are encoded in
  // an instance of AdjacentMatch 

#ifdef OLD_GAP
 void MatchAlgorithmGapped::generateMatches
  ( HitListVector& hitList, MatchAdder& addMatch ) 
  {

    //    subjectTable.matchSequence( querySeq, hitList, numRepeats_ );

    sort( hitList.begin(), hitList.end(), LessThanSubject() );

    //    bool foundMatch( false );

    if (hitList.size()==0) return;
    //    int origSize( size() );

    int numMatches( wordLength_ );

    HitListVector::iterator pFirstMatch( hitList.begin() );

    for ( int i(0) ; i < ((int)hitList.size()) - 1 ; ++i )
    {
      //                 cout << i << ": " 
      //   << (*this)[i].subjectNum << "|" 
      //   << (*this)[i].diff << "|" << (*this)[i].queryPos << " " 
      //   << (*this)[i+1].subjectNum << "|" 
      //   << (*this)[i+1].diff << "|" << (*this)[i+1].queryPos << " "
      //   << numMatches ; 

              //  && ( (*this)[i+1].queryPos > (*this)[i].queryPos )  )
            if (    ( hitList[i+1].subjectNum == hitList[i].subjectNum ) 
                 && ( hitList[i+1].diff - hitList[i].diff  <=maxGap_ )) 
      {
        //              cout << "+";
        if ( numMatches == 0 ) 
        { 
          pFirstMatch = &hitList[i]; 
          //                cout << "Start of run";
          numMatches = wordLength_; // %%%%
        }
        numMatches += wordLength_;
      } // ~if
      else if ( numMatches != 0 )
      {
        //      cout << "-";
          if ( numMatches >= minToProcess_ ) 
          { 
            //           cout << "storing match";


           sort( pFirstMatch, &hitList[i+1], LessThanQuery() );

           //        cout << "Storing match (" << numMatches 
           //         << " matches). Frequency array:\n";

           //            for ( HitListVectorGap::iterator z( pFirstMatch );
           //      z != &(*this)[i+1]; ++z )
           // cout<< z->subjectNum <<" "<< z->diff<< " "<< z->queryPos << "\n";

           //            foundMatch = true;

            addMatch ( pFirstMatch->subjectNum,
                      numMatches,
                      pFirstMatch->queryPos,
                      hitList[i].queryPos 
                      + wordLength_ - 1,
                      pFirstMatch->diff + pFirstMatch->queryPos,
                      hitList[i].diff + hitList[i].queryPos
                      + wordLength_ - 1 ); 
                    }
          numMatches = 0;
      } // ~else if
      //       cout << "\n";
    } // ~for

    if ( ( numMatches >= minToProcess_ ) && (!hitList.empty()) )
    {
        sort( pFirstMatch, hitList.end(), LessThanQuery() );
        // cout << "Storing match (" << numMatches 
        //  << " matches). Frequency array:\n";

        // for ( HitListVectorGap::iterator z( pFirstMatch );
        // z != end(); ++z )
        // cout<< z->subjectNum << " " << z->diff<< " " << z->queryPos<< "\n";

        //        foundMatch = true;

        addMatch ( pFirstMatch->subjectNum,
                    numMatches,
                    pFirstMatch->queryPos,
                    hitList.back().queryPos 
                    + wordLength_ - 1,
                    pFirstMatch->diff + pFirstMatch->queryPos,
                    hitList.back().diff + hitList.back().queryPos
                    + wordLength_ - 1 
                          ); 
    }

    return;// foundMatch; // ( size() != origSize );

  } // ~findMatch
#endif




// End of file QueryManagerSSAHA.cpp

---