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

// Description:

// Includes:

#include "HashTable.h"
#include "SequenceReader.h"
// #include "SequenceReaderCodon.h"
#include "GlobalDefinitions.h"
#include <fstream>
#include <algorithm>
#include <map>

AllocatorLocal<PositionInHitList> HashTable::defaultArrayAllocator;
AllocatorLocal<PositionInDatabase> HashTable::defaultHitListAllocator;

// ### Function Definitions ###



// Name:
// Arguments:
// TYPE  NAME  IN/OUT COMMENT
// Returns: TYPE COMMENT
  //  MatchPointer pM_;
  //  void fred
  //  ( WordSequence& seq, HitList& hitListFwd )
  //   { pM_=&HashTable::matchSequence;
  //     return (this->*pM_)(seq, hitListFwd); }

void HashTable::setNumRepeats( int numRepeats )
{
  if ( (numRepeats<0) || (numRepeats>stepLength_) )
    throw SSAHAException("Invalid value for numRepeats!!");
  numRepeats_=numRepeats;
  pMatchSequence_ = ( numRepeats==0 )
    ? &HashTable::matchSequenceStandard
    : &HashTable::matchSequenceRepeated;
}

  // Function Name: matchSequence
  // Arguments: WordSequence& (in), HitList& (out), HitList& (out)
  // Returns: void
  // This obtains the full list of hits for a sequence in both forward
  // and reverse directions. Proceeds as follows:
  // 1. The reverse complement of the sequence is formed.
  // 2. Any hits found in the forward or reverse direction are added to the
  // appropriate list.
  // 3. The sequence and reverse complement are left-shifted by 1 base
  // Steps 2 and 3 are repeated wordLength_ times.
  // NB This function will modify seq. If you want to keep it, make a copy
  // before calling this function.
  void HashTable::matchSequenceStandard
  ( WordSequence& seq, HitList& hitListFwd )
  {
    
    int numBasesInLast( seq.getNumBasesInLast() );

    for ( int i(0) ; i < wordLength_ ; ++i )
    {
      matchWord( seq,    hitListFwd, i );
      shiftSequence( seq, bitsPerSymbol_, wordLength_ );
      if ( i == numBasesInLast )
      {
        seq.pop_back();
      } // ~if
    } // ~for

  } // ~HashTable::matchSequence

  // Function Name: matchSequence
  // Arguments: WordSequence& (in), HitList& (out), HitList& (out), int (in)
  // Returns: void
  // This obtains the full list of hits for a sequence in both forward
  // and reverse directions and masks out tandem repeats.
  void HashTable::matchSequenceRepeated
  ( WordSequence& seq, 
    HitList& hitListFwd )
  {
    WordSequenceShifted seqShifted(seq, *this);
    //    screenRepeats( seqShifted, hitListFwd, numRepeats_ );

    Word                thisWord;

    int m;

    //  cout << "my size: " << size() << endl;

    // i cycles through each full word in the query sequence
    for ( int i(0) ; i < seqShifted.size() ; ++i )
    {


      //      cout << "doing  i:" << i << endl;
      thisWord = seqShifted[i];
      m = 0;

      // look through the next numRepeats_ words for duplicates
      for ( int j(i+1) ; 
            ( ( j < seqShifted.size() ) && ( j <= i + numRepeats_ ) ); 
            ++j )
      {
        if ( thisWord == seqShifted[j] )
        {
          m = j - i;
          //      cout << "Tandem repeat: " << i << "-" << j << "\n"; // %%%%
          break;
        } // ~if
      } // ~for j
      if ( m == 0 )
      {
        //      cout << "doing bog standard matching for:" << i << endl;
        matchWord( thisWord, hitListFwd, i );
      } // ~if
      else
      {
          // ... then we have found a tandem repeat of length m
          int r(1);

          // scan forward until we reach either the end of the
          // repeated region or the end of the sequence
          while (    ( seqShifted[i+(r*m)]==thisWord )
                  && ( i+(r*m) < seqShifted.size() ) ) ++r;

          //          cout << "Num repeats: " << r << endl;

          // any hits in a run of matching hits that exceed lastRun are
          // ignored, because in that case the size of the repeated
          // region in the subject sequence exceeds the size of the
          // repeated region in the query sequence
          int lastRun((r-1)*m);

          //      cout << "size of repeated run: " << lastRun << endl;

          while ( seqShifted[i+lastRun] == seqShifted[i+lastRun-m] ) lastRun++;
          lastRun--;

          //      cout << "adjusted size of repeated run: " << lastRun << endl;

          HitListRepeated hits;

          // as we proceed base by base along a region of tandem repeats
          // of motif length m, we encounter m distinct hash words, after
          // that, they repeat. Now get the hits for each of these words: 
          // passing in j tags each hit with its position in the repeat cycle
          for ( int j(0) ; j < m ; ++j ) 
          {
            matchWord( seqShifted[i+j], hits, j);
          } // ~for j

          // sort hits in order of RepeatedHit::subjectPos
          sort( hits.begin(), hits.end() );

          // lastHit =  previous hit in list, initialized to all zeroes
          RepeatedHit lastHit; 
          // firstHit = first hit of a matching run, initialized to all zeroes
          RepeatedHit firstHit;

          // thisRun = size of current run of matching hits
          int thisRun(0);

        
          for ( HitListRepeated::iterator thisHit( hits.begin() ); 
                thisHit != hits.end() ; ++thisHit )
          {
            //      cout << ": " << (*thisHit).subjectPos.sequence << " " 
            //      << (*thisHit).subjectPos.offset << " " 
            //      << (*thisHit).cyclePos ; 
            if (    (    (*thisHit).subjectPos.sequence 
                      == lastHit.subjectPos.sequence )
                 && (    (*thisHit).subjectPos.offset
                      == lastHit.subjectPos.offset + stepLength_ )
                 && (    (*thisHit).cyclePos
                      == (( lastHit.cyclePos + stepLength_ ) % m)   ) )
            {
              if ( thisRun == 0 ) 
              { 
                //      cout << " -s- ";
                // then a run if matching hits has started
                firstHit = lastHit; 
                thisRun  = wordLength_;
              } // ~if
              else 
              {
                //              cout << " -c- ";
                // we continue an existing run             
                thisRun += stepLength_;
              }
 

            } // ~if
            else 
            { 
              thisRun=0; 
              firstHit = *thisHit;
            }
            //      cout << " -n- ";
            if (thisRun <= lastRun ) 
            { 
              // only output hits if length of repeated region in subject 
              // is less than or equal to that of query
              //    cout << "added" << (*thisHit).subjectPos.offset 
              //   << "-" << i + firstHit.cyclePos + thisRun;
              hitListFwd.addHit( thisHit->subjectPos, 
                                   i + firstHit.cyclePos + thisRun ); 
            } else // cout << "ignored" << (*thisHit).subjectPos.offset 
                   //     << "-" << i + firstHit.cyclePos + thisRun;

       
            lastHit = *thisHit;
            //      cout << endl;
          } // ~for

          // carry on at the end of the repeated region in the query seq
          // i += (r-1)*m+1;
          // Actually want 1 less cos i is incremented anyway at the
          // end of the for loop - TC 5.7.1
          //          i += (r-1)*m;
          i += lastRun;
          // cout << "Carrying on at pos " << i+1 << "\n";
          // break out of the `for m' loop
          //     break;
       


      } // ~else

    } // ~for i



  } // ~HashTable::matchSequenceRepeated





void HashTable::countWords( SequenceAdapter& thisSeq )
{

   for ( int j(0) ; j < thisSeq.size() ; ++ j )
   {
     // only count words that have not been flagged
     pWordPositionInHitList_[(thisSeq[j]&(~gCursedWord))]
       += ((thisSeq[j]&gCursedWord)==(Word)0);
     //   pWordPositionInHitList_[thisSeq[j]]++;
   }
   
} // ~HashTable::countWords

void HashTable::hashWords
( SequenceAdapter& thisSeq, SequenceNumber seqNum )
{

        register Word              thisWord;
        register PositionInHitList currentPos;
      // NB We stop at the last but one element of the 
      // sequence (as the last isn't a full word)

   for ( int j(0) ; j < thisSeq.size() ; ++ j )
      {
        thisWord = thisSeq[j];
        // only hash words that have not been flagged
        if ((thisWord&gCursedWord)!=(Word)0) continue;
        currentPos 
        = pHitsFoundSoFar_[thisWord]
          +( ( thisWord == 0 ) 
             ? 0 : pWordPositionInHitList_[thisWord - 1]) ;
                           
        if ( currentPos != pWordPositionInHitList_[thisWord] )
        { // then place position in the hit list
          pHitListForAllWords_[currentPos].sequence 
          = seqNum; 
          pHitListForAllWords_[currentPos].offset   
          = j * stepLength_; 
          //      DEBUG_L2("list "<< printWord(thisWord,wordLength_) 
          //  << " "<< seqNum << " " << j*stepLength_ );
          pHitsFoundSoFar_[thisWord]++;
        } // ~if

      } // ~ for thisWord


} // ~HashTable::hashWords
  

// End of file HashTable.cpp

---