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

// Description:

// Includes:

#include "HashTablePacked.h"
#include "SequenceReader.h"
// #include "SequenceReaderCodon.h"
#include "GlobalDefinitions.h"
#include <fstream>
#include <algorithm>
#include <map>
#include <cmath> // for `pow'

AllocatorLocal<PositionInHitList> HashTablePacked::defaultArrayAllocator;
AllocatorLocal<PositionPacked> HashTablePacked::defaultHitListAllocator;

// ### Function Definitions ###

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

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

void HashTablePacked::setSubstituteThreshold( int numSubs ) 
{
  if ( (numSubs<0) || (numSubs>maxNumHits_) )
    throw SSAHAException("Invalid value for substituteThreshold!!");
  substituteThreshold_=numSubs;
  pMatchWord_ = ( numSubs==0 )
    ? &HashTablePacked::matchWordStandard
    : &HashTablePacked::matchWordSubstitute;
} // ~HashTablePacked::setSubstituteThreshold( int numSubs ) 

void HashTablePacked::matchWordSubstitute
( Word w, PackedHitStore& hitList, int offset )
{
  if ((w&gCursedWord)!=(Word)0) return;

  HashTableView<PositionPacked,HashTablePacked>::matchWord
    ( w, hitList, offset );
  assert((w&gCursedWord)==0);
  if (size(w)>substituteThreshold_) return;
  vector<Word> neighbours;
  (*pGenerateSubstitutes_)( w, neighbours, wordLength_);
  for( vector<Word>::iterator i(neighbours.begin());
       i!=neighbours.end(); i++)
  {
    HashTableView<PositionPacked,HashTablePacked>::matchWord
      ( *i, hitList, offset );
  } // ~for i
} // void HashTablePacked::matchWordSubstitute

// Given a DNA word, generates all single base purine-purine (A-G) and
// pyrimidine-pyrimidine (C-T) substitutes. 
// A = 00, A XOR 10 = 10 = G
// C = 01, C XOR 10 = 11 = T
// G = 10, G XOR 10 = 00 = A
// T = 11, T XOR 10 = 01 = C
void generateSubstitutesDNA
(Word w, vector<Word>& subs, int wordLength)
{
  Word mask(2);
  for (int i(0) ; i < wordLength ; i++, mask<<=gBaseBits)
  {
    subs.push_back(w^mask);
  } // ~for
} // ~generateSubstitutesDNA

// Given a protein word, generates all single residue substitutes having a
// positive BLOSUM62 score.
void generateSubstitutesProtein
(Word w, vector<Word>& subs, int wordLength)
{
  Word mask((1<<gResidueBits)-1), thisCodon, thisWord;
  for (int i(0) ; i < wordLength ; i++, mask<<=gResidueBits)
  {
    thisCodon=((w&mask)>>(i*gResidueBits));
    thisWord= (w&(~mask));
    //    cout << ~mask << " " << printResidue( thisWord, wordLength )
    // << " " << printResidue( thisCodon, wordLength) << endl;
    for (int j(subStarts[thisCodon]); j != subStarts[thisCodon+1]; j++ )
    {
      subs.push_back( thisWord | (subVals[j]<<(i*gResidueBits)) );
      //  cout << "generating sub: " << printResidue(subs.back(),wordLength)
      //  << endl;
    } // ~for j
  } // ~for i
} // ~generateSubstitutesProtein





void HashTablePacked::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]]++;
   }   
   //   cout << endl;
} // ~HashTable::countWords

void HashTablePacked::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)
        {

          currentPos 
            = pHitsFoundSoFar_[thisWord]
            +( ( thisWord == 0 ) 
               ? 0 : pWordPositionInHitList_[thisWord - 1]) ;
                           
          if ( currentPos != pWordPositionInHitList_[thisWord] )
          { // then place position in the hit list
            pHitListForAllWords_[currentPos] = wordNum_;
            pHitsFoundSoFar_[thisWord]++;
            // next line moved: wordNum still needs incrementing even if
            // word is flagged
          //      wordNum_++;
          } // ~if
          else assert(1==0);
        } // ~if

        wordNum_++;

      } // ~ for thisWord

   seqStarts_.push_back(wordNum_);

} // ~HashTable::hashWords

void HashTablePacked::matchSequenceStandard
( WordSequence& seq, HitList& hitListFwd )
{

  PackedHitStore packedHits;
  int numBasesInLast( seq.getNumBasesInLast() ), baseOffset;

  for ( int i(0) ; i < wordLength_ ; ++i )
  {
    //  matchWord( seq,    packedHits, i );
    if (seq.size()!=0)
    {
      baseOffset=i;
      WordSequence::const_iterator last(&seq.back());
      for ( WordSequence::const_iterator thisWord(seq.begin());
            thisWord != last ; ++thisWord )
      {
        int oldSize(packedHits.size()); // %%%%%%
        matchWordDeluxe( *thisWord, packedHits, baseOffset );
        //      cout << printResidue(*thisWord, wordLength_) << " "
        //  << packedHits.size()-oldSize;
        //      for (int fk(oldSize);fk!=packedHits.size();fk++) cout << " " << packedHits[fk].first;
        //      cout << endl;
        baseOffset += wordLength_;
      } // ~for
    } // ~if
    shiftSequence( seq, bitsPerSymbol_, wordLength_ );
    if ( i == numBasesInLast )
    {
      seq.pop_back();
    } // ~if
  } // ~for

  convertHits(packedHits,hitListFwd);
  
} // ~HashTablePacked::matchSequenceStandard


void HashTablePacked::convertHits
( PackedHitStore& packedHits, HitList& hitListFwd )
{

  //  sort(packedHits.begin(),packedHits.end());
  sorter_(packedHits);


  vector<HitPacked>::iterator i(packedHits.begin());
  vector<SeqStartPos>::iterator ub, j(seqStarts_.begin());
  HitListVector::size_type sortStart;
  //  hitListFwd.reserve(hitListFwd.size()+packedHits.size());

  // NB seqStarts_ must not be empty else call to back() segfaults

  while 
  (    (j!=static_cast<vector<SeqStartPos>::iterator>(&seqStarts_.back()))
    && (i!=packedHits.end()) )
  {
    ub = upper_bound(j,seqStarts_.end(),i->first);
    j=ub; j--;
    sortStart=hitListFwd.size();
    while ((i!=packedHits.end())&&(i->first<*ub))
    {
      //      cout << i->first << " ... " 
      //   << ub-seqStarts_.begin() << " " 
      //   << stepLength_*(i->first - *j) 
      //   << endl;
      hitListFwd.addHit( ub-seqStarts_.begin(), 
                         stepLength_*(i->first - *j),
                         i->second );
      i++;
      
    } // ~while
    sort
    ( static_cast<HitList::iterator>(&hitListFwd[sortStart]),
      hitListFwd.end(),
      LessThanDiff() );
    //    sort(hitListFwd.begin()+sortStart,hitListFwd.end(),LessThanDiff() );
    j=ub;  
  } // ~while

} // ~HashTablePacked::convertHits

void HashTablePacked::matchSequenceRepeated
( WordSequence& seq, HitList& hitListFwd )
{
    WordSequenceShifted seqShifted(seq, *this);
    //    screenRepeats( seqShifted, hitListFwd, numRepeats_ );

    PackedHitStore nonRepeatedHits;

    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;
        matchWordDeluxe( thisWord, nonRepeatedHits, 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;
          PackedHitStore packedHits;
          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], packedHits, j);
          } // ~for j

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

          vector<HitPacked>::iterator pThisHit(packedHits.begin());
          vector<SeqStartPos>::iterator ub, j(seqStarts_.begin());
          //  hitListFwd.reserve(hitListFwd.size()+packedHits.size());
          
          // NB seqStarts_ must not be empty else call to back() segfaults

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

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

          while 
          (    ( j!=static_cast<vector<SeqStartPos>::iterator>
                      ( &seqStarts_.back() ) ) 
            && ( pThisHit!=packedHits.end() ) )
          {
            ub = upper_bound(j,seqStarts_.end(),pThisHit->first);
            j=ub; j--;
            lastHit=*pThisHit; 
            // ensures if condition always false first time through while
            while ((pThisHit!=packedHits.end())&&(pThisHit->first<*ub))
            {
              if (    (    pThisHit->first
                           == lastHit.first+1 )
                      && (    pThisHit->second
                              == (( lastHit.second + stepLength_ ) % m)   ) )
              {
                if ( thisRun == 0 ) 
                { 
                  //    cout << " -s- ";
                  // then a run of matching hits has started
                  firstHit = lastHit; 
                  thisRun  = wordLength_;
                } // ~if
                else 
                {
                  //            cout << " -c- ";
                  // we continue an existing run             
                  thisRun += stepLength_;
                } // ~else
 

              } // ~if
              else 
              { 
                thisRun=0; 
                firstHit = *pThisHit;
              }
              //            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 ); 

                nonRepeatedHits.push_back
                  ( HitPacked( pThisHit->first, 
                               i + firstHit.second + thisRun ) ); 
              } // ~if
              else 
              {
                lastHit = *pThisHit;
              } // ~else
              pThisHit++;
            } // ~while
            j=ub;  
          } // ~while
          i += lastRun;
      } // ~else
    } // ~for i
    convertHits( nonRepeatedHits, hitListFwd );
    
} // ~HashTablePacked::matchSequenceRepeated




void HashTablePacked::loadHitList( unsigned long size )
{
  
  string startFileName(name_+(string)".start");

  ifstream startFile(startFileName.c_str());
  if (startFile.fail())
  {
    monitoringStream_ << "Could not open " << startFileName << endl;
    throw SSAHAException("Could not open .start file");
  } // ~if

  startFile.seekg(0,ios::end);
  long startFileSize = startFile.tellg();

  if ( getNumSequences() + 1 != startFileSize/sizeof(SeqStartPos))
  {
    seqStarts_.resize(startFileSize/sizeof(SeqStartPos));
    monitoringStream_ 
      << "Info: expecting " << seqStarts_.size()-1 << " sequences in file"
      << endl;
  } // ~if
  else seqStarts_.resize(getNumSequences()+1);

  loadFromFile(  startFileName,
                 (char*) &seqStarts_[0],
                 startFileSize  );

  HashTableView<PositionPacked,HashTablePacked>::loadHitList(size);

}



void HashTablePacked::saveHitList( void )
{

  //  assert(seqStarts_.size()==getNumSequences()+1);
  // assertion removed, not true when part of a HashTableTranslated

  saveToFile(  name_+(string)".start",
               (char*) &seqStarts_[0],
               (seqStarts_.size())*sizeof(SeqStartPos)  );

  HashTableView<PositionPacked,HashTablePacked>::saveHitList();
} // ~HashTablePacked::saveHitList( void )


// RadixSorter member function definitions
RadixSorter::RadixSorter
( const unsigned int digits, const unsigned int bits ) :
digits_(digits), bits_(bits), 
base_( (unsigned long)pow( (double)2, (int)bits ) ), 
mask_( base_ - 1 ),
counts_(vector< vector<CountInt> >(digits, vector<CountInt>(base_) ) ),
places_(base_)
{} // RadixSorter::RadixSorter



void RadixSorter::operator()( vector<HitPacked>& v )
{

  CountDigits(v);

  if (digits_&1) 
  { 
    v1_=v;
    source_ = &v1_; target_ = &v;
  } // ~if
  else 
  { 
    v1_.assign(v.size(),zeroHit);
    source_ = &v; target_ = &v1_;
  } // ~else

  vector<HitPacked>* temp;

  for ( unsigned int i(0) ; i < digits_ ; ++i ) 
  {
    SortByDigit( i );
    temp = source_;
    source_=target_;
    target_=temp;
  } // ~for
} // ~void RadixSorter::operator()( vector<HitPacked>& v )


void RadixSorter::CountDigits( const vector<HitPacked>& v )
{
  vector< vector<CountInt> >::iterator pv (counts_.begin()); 

  for (; pv != counts_.end() ; ++pv ) pv->assign(base_,0);

  PositionPacked temp;

  for ( vector<HitPacked>::const_iterator 
          j(v.begin()); 
          j != v.end() ; ++j )
  {
    temp = j->first;
    for ( pv = counts_.begin() ; pv != counts_.end() ; ++pv )
    {
      (*pv)[ temp & mask_ ]++; 
      temp >>= bits_;
    } // ~for pv

  } // ~for j

} // ~RadixSorter::CountDigits( const vector<Uint>& v )

void RadixSorter::SortByDigit( PositionPacked digit )
{
  unsigned int shift( digit * bits_ );
  CountInt pos(0);
  vector<CountInt>::iterator count( counts_[digit].begin() ); 

  for ( vector< vector<HitPacked>::iterator >::iterator 
          place( places_.begin() );
        place != places_.end(); 
        ++place, ++count )
  {
    //    *place = &((*target_)[pos]);
    *place = static_cast<vector<HitPacked>::iterator>(&(*target_)[pos]);
    //    *place = target_->begin()+pos;
    pos += *count;
  } // ~for

  for ( vector<HitPacked>::iterator i( (source_)->begin() ) ; 
        i!= (source_)->end() ; ++i )
  {
    *places_[ ( i->first >> shift) & mask_ ]++ = *i;
  } // ~for
} // ~void RadixSorter::SortByDigit( PositionPacked digit )


// End of file HashTablePacked.cpp

---