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

// Description:

// Includes:

#include "HashTableTranslated.h"
#include "SequenceReader.h"
#include "SequenceEncoder.h"

// ### Function Definitions ###

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




// HashTableComponent member function definitions

HashTableComponent::HashTableComponent( ostream& monitoringStream, string name,
                   Allocator<PositionPacked>& hitListAllocator,
                   Allocator<PositionInHitList>& arrayAllocator 
                   ):
  queryFrame_(0),
  HashTablePacked( monitoringStream, name, 
                   hitListAllocator, arrayAllocator  ) 
{
  hitListFormat_ = g32BitPackedProtein;
  bitsPerSymbol_ = gResidueBits;
  pGenerateSubstitutes_ = &generateSubstitutesProtein;
  monitoringStream_ << "constructing HashTableComponent" << endl;
  pNameReader_ = new NameReader; // dummy, component contains no name info
} // ~constructor



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

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

  vector<HitPacked>::iterator i(packedHits.begin());
  vector<SeqStartPos>::iterator j(seqStarts_.begin());
  HitListVector::size_type sortStart;

  SequenceNumber adjustedSeqNum;
  int subjectFrame;

  //  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()) )
  {
    vector<SeqStartPos>::iterator ub(upper_bound(j,seqStarts_.end(),i->first));
    j=ub; j--;

    sortStart=hitListFwd.size();
    while ((i!=packedHits.end())&&(i->first<*ub))
    {
      //       cout << "Hit: " << i->first << " ... " 
      //   << ub-seqStarts_.begin() << " " 
      //  << stepLength_*(i->first - *j) << " "
      // << i->second      << endl;

      adjustedSeqNum 
        = ub-seqStarts_.begin() - 1; // 0,1,2 ... etc.
      subjectFrame = adjustedSeqNum % gNumReadingFrames;
      adjustedSeqNum /= gNumReadingFrames;
      adjustedSeqNum++;


      //    cout 
      //        << "aHit: "
      //<< adjustedSeqNum << " " 
      //    << (gNumReadingFrames*stepLength_*(i->first - *j))+subjectFrame << " "
      //    << (gNumReadingFrames*i->second)+queryFrame_ << endl; 
        hitListFwd.addHit
          ( adjustedSeqNum,
            (gNumReadingFrames*stepLength_*(i->first - *j))+subjectFrame,
            (gNumReadingFrames*i->second)+queryFrame_); 
        //      hitListFwd.addHit( ub-seqStarts_.begin(), 
        //       stepLength_*(i->first - *j),
        //       i->second );
      i++;
      
    } // ~while i
    sort
    ( static_cast<HitList::iterator>(&hitListFwd[sortStart]),
      hitListFwd.end(),
      LessThanDiff() );
    j=ub;  
  } // ~while j

} // ~HashTableComponent::convertHits

// HashTableProteinPacked member function definitions


HashTablePackedProtein::HashTablePackedProtein
( ostream& monitoringStream, string name,
  Allocator<PositionPacked>& hitListAllocator,
  Allocator<PositionInHitList>& arrayAllocator 
  ):
  queryFrame_(0),
  queryMult_(1),
  codonEncoder_(5),
  pMatchSequence_( &HashTablePackedProtein::matchSequenceProtein ),
  HashTablePacked( monitoringStream, name, 
                     hitListAllocator, arrayAllocator  ) 
{
    hitListFormat_ = g32BitPackedProtein;
    bitsPerSymbol_ = gResidueBits;
    pGenerateSubstitutes_ = &generateSubstitutesProtein;
    monitoringStream_ << "constructing HashTablePackedProtein" << endl;
} // constructor





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


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

  vector<HitPacked>::iterator i(packedHits.begin());
  vector<SeqStartPos>::iterator j(seqStarts_.begin());
  HitListVector::size_type sortStart;

  //SequenceNumber adjustedSeqNum;
  //int subjectFrame;

  //  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()) )
  {
    vector<SeqStartPos>::iterator ub(upper_bound(j,seqStarts_.end(),i->first));
    j=ub; j--;

    sortStart=hitListFwd.size();
    while ((i!=packedHits.end())&&(i->first<*ub))
    {
      //            cout << "Hit: " << i->first << " ... " 
      // << ub-seqStarts_.begin() << " " 
      // << stepLength_*(i->first - *j) << " "
      // << i->second      << endl;

            //    cout 
            //  << "aHit: "
            //  << ub-seqStarts_.begin() << " "
            //  << queryMult_*stepLength_*(i->first - *j) << " "
            //<< (queryMult_*i->second)+queryFrame_ << endl; 
        hitListFwd.addHit
          ( ub-seqStarts_.begin(), 
            (queryMult_*stepLength_*(i->first - *j)),
            (queryMult_*i->second)+queryFrame_); 
        //      hitListFwd.addHit( ub-seqStarts_.begin(), 
        //       stepLength_*(i->first - *j),
        //       i->second );
      i++;
      
    } // ~while i
    sort
    ( static_cast<HitList::iterator>(&hitListFwd[sortStart]),
      hitListFwd.end(),
      LessThanDiff() );
    j=ub;  
  } // ~while j


} // ~HashTablePackedProtein::convertHits

void HashTablePackedProtein::matchSequenceProtein
( WordSequence& seq, HitList& hitListFwd )
{
  HashTablePacked::matchSequence( seq, hitListFwd );
} // ~HashTablePackedProtein::matchSequenceProtein


void HashTablePackedProtein::matchSequenceTranslatedDNA
( WordSequence& seq, HitList& hitListFwd )
{
  codonEncoder_.setWordLength( wordLength_ );

  WordSequence translatedSeq;
  for (queryFrame_ = 0; queryFrame_ < gNumReadingFrames; queryFrame_++)
  {
    codons_.clear();
    codonize ( seq, codons_, queryFrame_ );
    translatedSeq.clear();
    codonEncoder_.linkSeq( translatedSeq );
    codonEncoder_.encode( codons_ );
    codonEncoder_.unlinkSeq();
    HashTablePacked::matchSequence( translatedSeq, hitListFwd );
  } // ~for

} // ~HashTablePackedProtein::matchSequenceTranslatedDNA






// HashTableTranslated member function definitions

HashTableTranslated::HashTableTranslated
( ostream& monitoringStream, string name,
  Allocator<PositionPacked>& hitListAllocator,
  Allocator<PositionInHitList>& arrayAllocator
  ) :
  hashFwd_( monitoringStream, name+(string)"_fwd", 
            hitListAllocator, arrayAllocator ),
  hashRev_( monitoringStream, name+(string)"_rev", 
            hitListAllocator, arrayAllocator ),
  pHash_(&hashFwd_),
  codonEncoder_(5),
  pMatchSequence_( &HashTableTranslated::matchSequenceProtein ),
  HashTableGeneric( monitoringStream, name, arrayAllocator ) 
{
  bitsPerSymbol_=gResidueBits;
  hitListFormat_ = gTranslated;
  monitoringStream_ << "constructing HashTableTranslated" << endl;
} // ~constructor

void HashTableTranslated::setupPointerArray( void )
{

  hashFwd_.numDifferentWords_= 1 << (wordLength_*gResidueBits);
  hashFwd_.bitsPerSymbol_=gResidueBits;
  hashFwd_.sourceData_=gDNAData;
  hashFwd_.wordLength_=wordLength_;
  hashFwd_.stepLength_=stepLength_;
  hashFwd_.setupPointerArray();

  hashRev_.numDifferentWords_= 1 << (wordLength_*gResidueBits);
  hashRev_.bitsPerSymbol_=gResidueBits;
  hashRev_.sourceData_=gDNAData;
  hashRev_.wordLength_=wordLength_;
  hashRev_.stepLength_=stepLength_;
  hashRev_.setupPointerArray();
}


void HashTableTranslated::computePointerArray( void )
{
  hashFwd_.computePointerArray();
  hashRev_.computePointerArray();
}



void HashTableTranslated::setupHitList( void )
{
  hashFwd_.setupHitList();
  hashRev_.setupHitList();
}

void HashTableTranslated::cleanupTempData( void )
{
  hashFwd_.cleanupTempData();
  hashRev_.cleanupTempData();
  hashFwd_.isInitialized_=true;
  hashRev_.isInitialized_=true;
}


void HashTableTranslated::hashWords
( SequenceAdapter& thisSeq, SequenceNumber seqNum ) 
{
  assert(1==0);
  hashFwd_.hashWords( thisSeq, seqNum );
  hashRev_.hashWords( thisSeq, seqNum );
}

int  HashTableTranslated::getMaxNumHits() const
{ 
  assert( hashFwd_.getMaxNumHits() == hashRev_.getMaxNumHits() );
  return hashFwd_.getMaxNumHits();
} // ~getMaxNumHits() const 
void HashTableTranslated::setMaxNumHits( int mnh )
{ 
  hashFwd_.setMaxNumHits( mnh );
  hashRev_.setMaxNumHits( mnh );
  assert( hashFwd_.getMaxNumHits() == mnh );
  assert( hashFwd_.getMaxNumHits() == mnh );
} // ~setMaxNumHits( int mnh ) 



int HashTableTranslated::countWordsAndGetNames
( SequenceReader& sequenceReader, SequenceAdapter* seq )
{

    sequenceReader.rewind();

    //    CodonList codons; // make class member %%%%%
    //  SequenceEncoderCodon codonEncoder(5); // make class member %%%%%
    //    cout << "WordLength :" << wordLength_ << endl;
    NameReaderLocal* pReaderLocal = new NameReaderLocal;
    pNameReader_ = pReaderLocal;

    codonEncoder_.setWordLength(wordLength_);

    // Change encoding mode
    SequenceReaderModeFlagReplace mode('X');
    assert(ttCodon['X']==ttProtein['X']);
    assert(ttCodon['X']!=nv);
    codonEncoder_.changeMode( &mode );


    //    WordSequence thisSeq, revSeq, translatedSeq;
    WordSequence thisSeq, translatedSeq;
    seq->link( translatedSeq );
    int numSeqs(0);

    //    sequenceNames_.push_back(new string);


    while( sequenceReader.getNextSequence
           ( thisSeq, eDNAWordSizeForHashing ) != -1 ) 
    {
      numSeqs++;

      //      sequenceReader.getLastSequenceName( *(sequenceNames_.back()) );
      sequenceReader.getLastSequenceName( pReaderLocal->lastName() );
      //      sequenceNames_.push_back( new string );

      //    hashRev_.countWords(*seq);

      for ( int i(0) ; i < gNumReadingFrames ; i++ )
      {
        codons_.clear();
        codonizeAndFlag ( thisSeq, codons_, i );
        translatedSeq.clear();
        codonEncoder_.linkSeq( translatedSeq );
        codonEncoder_.encode( codons_ );
        codonEncoder_.unlinkSeq();
        hashFwd_.countWords(*seq);
      } // ~for i
      
      //      revSeq.clear();
      //    reverseComplement( thisSeq, revSeq, eDNAWordSizeForHashing);

      for ( int i(0) ; i < gNumReadingFrames ; i++ )
      {
        codons_.clear();
        codonizeAndFlagReverse ( thisSeq, codons_, i );
        translatedSeq.clear();
        codonEncoder_.linkSeq( translatedSeq );
        codonEncoder_.encode( codons_ );
        codonEncoder_.unlinkSeq();
        hashRev_.countWords(*seq);
      } // ~for i


    } // ~while

    //    pReaderLocal->pop_back();
    //   pNameReader_ = pReaderLocal;

    //    delete sequenceNames_.back();
    //   sequenceNames_.pop_back();
    //    wordSeqs.pop_back();

    return numSeqs;
} // ~HashTableTranslated::countWordsAndGetNames

void HashTableTranslated::hashAllWords
( SequenceReader& sequenceReader, SequenceAdapter* seq, int numSeqs )
{

  //  cout << "HTT::hashAllWords" << endl;
  int numWords(0);
  sequenceReader.rewind();
  //  WordSequence thisSeq, revSeq, translatedSeq;
  WordSequence thisSeq, translatedSeq;
  seq->link( translatedSeq );

  // NB sequences are numbered 1...n not 0...n-1
  for ( unsigned int i(1); i <= numSeqs ; i++ )
  {
    //    cout << "hashing sequence " << i << endl;
    if( sequenceReader.getNextSequence( thisSeq, eDNAWordSizeForHashing) == -1 ) 
    {
      throw SSAHAException
        ("Sequence source data changed during hash table creation!");
    }

    // fill in the word count
    numWords = (int) thisSeq.size();
    pSequenceSizes_[i-1] = ( numWords > 0 ) 
      ? ( (numWords-1) * eDNAWordSizeForHashing ) + thisSeq.getNumBasesInLast()
      : 0;

    //    cout << "Fwd: " << printWord(thisSeq, eDNAWordSizeForHashing ) << endl;

    for ( int j(0) ; j < gNumReadingFrames ; j++ )
    {
      codons_.clear();
      codonizeAndFlag ( thisSeq, codons_, j );
      translatedSeq.clear();
      codonEncoder_.linkSeq( translatedSeq );
      codonEncoder_.encode( codons_ );
      codonEncoder_.unlinkSeq();

      //
      //            cout << "T" << j << ": ";
      //      for (WordSequence::iterator b(translatedSeq.begin());
      //   b!=translatedSeq.end();b++)
      //        cout << printResidue(*b, wordLength_ ) <<"-";
      //  cout << endl;


      hashFwd_.hashWords(*seq,99999);
    } // ~for
      
    //    revSeq.clear();
    //    reverseComplement( thisSeq, revSeq, eDNAWordSizeForHashing);
    
    //  cout << "Rev: " << printWord(thisSeq, eDNAWordSizeForHashing ) << endl;

    for ( int j(0) ; j < gNumReadingFrames ; j++ )
    {
      codons_.clear();
      codonizeAndFlagReverse ( thisSeq, codons_, j );
      translatedSeq.clear();
      codonEncoder_.linkSeq( translatedSeq );
      codonEncoder_.encode( codons_ );
      codonEncoder_.unlinkSeq();
      //           cout << "T" << j << ": ";
      //  for (WordSequence::iterator b(translatedSeq.begin());
      //    b!=translatedSeq.end();b++)
      //        cout << printResidue(*b, wordLength_ ) <<"-";
      //  cout << endl;

      hashRev_.hashWords(*seq,99999); // seqNum not needed for HTPacked::hash
    }

    // do the hashing
    //   seq->link(thisSeq);
    //   hashWords( *seq,i );
    //   cout << "hashallwords: " << wordLength_ << getWordLength() << endl;

  } // ~for thisSeq

  SequenceReaderModeIgnore mode;
  codonEncoder_.changeMode( &mode );


} // ~HashTableGeneric::hashAllWords 



void HashTableTranslated::countWords( SequenceAdapter& thisSeq ) 
{
  assert(1==0);
  //  hashFwd_.countWords( thisSeq );
  //  hashRev_.countWords( thisSeq );
}


void HashTableTranslated::setNumRepeats( int nr) 
{
  hashFwd_.setNumRepeats( nr );
  hashRev_.setNumRepeats( nr );
}

void HashTableTranslated::setSubstituteThreshold( int ns) 
{
  hashFwd_.setSubstituteThreshold( ns );
  hashRev_.setSubstituteThreshold( ns );
}

char* HashTableTranslated::getHitListStart( void ) const 
{
  assert(1==0);
  return NULL;
}
int HashTableTranslated::getHitTypeSize( void ) const 
{
  assert(1==0);
  return NULL;
}
void HashTableTranslated::allocateHitList( unsigned long size ) 
{
  assert(1==0);
}
void HashTableTranslated::loadHitList( unsigned long size ) 
{
  assert(1==0);
}
void HashTableTranslated::saveHitList( void ) 
{
  //  assert(1==0);
  monitoringStream_ << "No hit list to save for HashTableTranslated" << endl;
}

void HashTableTranslated::loadHashTable
( SourceReaderIndex* pSourceReader )
{

  if (pSourceReader==NULL) 
  {
    pNameReader_ = new NameReaderLocal;
  } // ~if
  else
  {
    pNameReader_ = new NameReaderIndex(*pSourceReader);
  } // ~else

  // These have their own dummy NameReaders
  hashFwd_.loadHashTable();
  hashRev_.loadHashTable();
  loadSequenceNames();

  //  unsigned long numSeqs = sequenceNames_.size();
  unsigned long numSeqs = pNameReader_->size();

  assert(numSeqs!=0);

  monitoringStream_ << "Allocating memory for pSequenceSizes_: "
                    << numSeqs << " sequences, "
                    << numSeqs*sizeof(SequenceOffset)
                    << " bytes total ...\n";

  pSequenceSizes_ = new SequenceOffset [ numSeqs  ];
      
  if (!pSequenceSizes_) 
  {
    throw SSAHAException("Memory allocation failed!");
  } // ~if

  loadFromFile(name_+(string)".size", (char*)pSequenceSizes_, 
               numSeqs * sizeof( SequenceOffset ),
               monitoringStream_ );

  codonEncoder_.setWordLength( wordLength_ );

  isInitialized_=true;
}


void HashTableTranslated::saveHashTable( void )
{
  hashFwd_.saveHashTable();
  hashRev_.saveHashTable();
  HashTableGeneric::saveHashTable();
  //  saveSequenceNames();
}

unsigned long HashTableTranslated::getTotalNumWords( void ) const
{
  return (hashFwd_.getTotalNumWords()+hashRev_.getTotalNumWords());
} // ~HashTableTranslated::getTotalNumWords( void ) const

void HashTableTranslated::printHashStats( void )
{
  cout << "\n\n ** Hash stats output for forward subtable: ** \n\n";
  hashFwd_.printHashStats();
  cout << "\n\n ** Hash stats output for reverse subtable: ** \n\n";
  hashRev_.printHashStats();
} // ~HashTableTranslated::printHashStats( void )



void HashTableTranslated::matchSequenceProtein
( WordSequence& seq, HitList& hitListFwd )
{
  //cout << "Trans::matchSequenceProtein, wordLength=" << wordLength_ << endl;
  //cout << printResidue( seq, wordLength_ ) << endl;
  pHash_->setQueryFrame(0);
  pHash_->matchSequence( seq, hitListFwd );
} // ~HashTableTranslated::matchSequenceProtein

void HashTableTranslated::matchSequenceTranslatedDNA
( WordSequence& seq, HitList& hitListFwd )
{
  //cout << "Trans::matchSequenceTranslatedDNA, wordLength=" << wordLength_ 
  //<< endl;
  codonEncoder_.setWordLength(wordLength_);
  WordSequence translatedSeq;
  for (int readingFrame(0); readingFrame < gNumReadingFrames; readingFrame++)
  {
    pHash_->setQueryFrame(readingFrame);
    codons_.clear();
    codonize ( seq, codons_, readingFrame );
    translatedSeq.clear();
    codonEncoder_.linkSeq( translatedSeq );
    codonEncoder_.encode( codons_ );
    codonEncoder_.unlinkSeq();
    //      cout << "T" << readingFrame << ": ";
    //  for (WordSequence::iterator b(translatedSeq.begin());
    //   b!=translatedSeq.end();b++)
    //       cout << printResidue(*b, wordLength_ ) <<"-";
    //   cout << endl;

    pHash_->matchSequence( translatedSeq, hitListFwd );
  } // ~for
} // ~HashTableTranslated::matchSequenceTranslatedDNA



// End of file HashTableTranslated.cpp

---