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

// Description:

// #pragma implementation
// Includes:

#include "SequenceReaderMulti.h"

// ### Function Definitions ###

// Name:
// Arguments:
// TYPE  NAME  IN/OUT COMMENT
// Returns: TYPE COMMENT
SequenceReaderMulti::SequenceReaderMulti( ostream& monitoringStream )
  : SequenceReader( monitoringStream ), thisReader_( allReaders_.begin() ),
    isFirstSeq_(true)
{
  monitoringStream_ << "constructing SequenceReaderMulti" << endl;
}

SequenceReaderMulti::~SequenceReaderMulti()
{
  monitoringStream_ << "destructing SequenceReaderMulti" << endl;
  for ( vector<SeqReaderInfo>::iterator i( allReaders_.begin() );
        i != allReaders_.end() ; i++ ) delete i->ptr_;
} // destructor

SequenceReaderMulti::SequenceReaderMulti( const SequenceReaderMulti& rhs) 
  : SequenceReader( rhs.monitoringStream_ ),
    allReaders_( rhs.allReaders_ ),
    thisReader_( allReaders_.begin() ),
    isFirstSeq_( rhs.isFirstSeq_ ),
    bitsPerSymbol_( rhs.bitsPerSymbol_ ),
    sourceDataType_( rhs.sourceDataType_ )
{
  monitoringStream_ << "copy constructing SequenceReaderMulti\n";
  for ( vector<SeqReaderInfo>::const_iterator i(rhs.allReaders_.begin());
        i!=rhs.allReaders_.end(); i++)
  {
    allReaders_.push_back( SeqReaderInfo() );
    allReaders_.back().ptr_=i->ptr_->clone();
  } // ~for

} // ~constructor


  
//template< class T >
//void SequenceReaderMulti::addSequenceReader( const T& seq )
//{
//allReaders_.push_back
//( std::auto_ptr<SequenceReader>( new T( seq ) ) );
//} // ~addSequenceReader( const SequenceReader& pSeq )


void SequenceReaderMulti::rewind( void )
{

    for ( vector<SeqReaderInfo>::iterator i 
          = allReaders_.begin(); 
          i != allReaders_.end();
          i++ )
    {
      (i->ptr_)->rewind();       
    } // ~for i
    thisReader_ = allReaders_.begin();
    lastSequenceNumber_ = 0;
} // ~SequenceReaderMulti::rewind( void )

  // Function Name: addReader
  // Arguments: const T&
  // A copy of seq is made (which is why SequenceReader and its subclasses
  // need to have a copy constructor) and a pointer to it is placed in 
  // allReaders. NB Type safety is ensured because if T is anything other than 
  // a subclass of SequenceReader then trying to 'push_back' a pointer to it
  // onto allReaders (a vector of SequenceReader pointers) will cause
  // the compiler to complain.
  //  template< class T >
  void SequenceReaderMulti::addReader( SequenceReader& seq )
  {
    if (isFirstSeq_)
    {
      isFirstSeq_ = false;
      bitsPerSymbol_ = seq.getBitsPerSymbol();
      sourceDataType_ = seq.getSourceDataType();
    } // ~if
    else
    {
      assert(bitsPerSymbol_ == seq.getBitsPerSymbol());
      assert(sourceDataType_ == seq.getSourceDataType());
    } // ~else
    //    assert(1==0);
    //    allReaders_.push_back( SeqReaderInfo( seq.clone() ));
    allReaders_.push_back( SeqReaderInfo() );
    allReaders_.back().ptr_=seq.clone();
    thisReader_
      = static_cast<vector<SeqReaderInfo>::iterator>(&allReaders_.back());
  } // ~addReader( SequenceReader& seq )

  // Function Name: addReader
  // Arguments: SequenceReader*
  // pSeq is placed in allReaders_. NB no copy of pSeq is made: use by
  // creating a new SequenceReader and calling this. The SeqReaderMulti
  // takes ownership of *pSeq and is responsible for its destruction
  void SequenceReaderMulti::addReader( SequenceReader* pSeq )
  {
    if (isFirstSeq_)
    {
      isFirstSeq_ = false;
      bitsPerSymbol_ = pSeq->getBitsPerSymbol();
      sourceDataType_ = pSeq->getSourceDataType();
    } // ~if
    else
    {
      assert(bitsPerSymbol_ == pSeq->getBitsPerSymbol());
      assert(sourceDataType_ == pSeq->getSourceDataType());
    } // ~else

    allReaders_.push_back( SeqReaderInfo() );
    allReaders_.back().ptr_=pSeq;
    thisReader_
      = static_cast<vector<SeqReaderInfo>::iterator>(&allReaders_.back());
  } // ~addReader( SequenceReader* pSeq )

// Function Name: findSequence
// Arguments: SequenceNumber (in)
// Returns:   void
// Winds the input file stream to the start of sequence number seqNum. 
// Returns false if seqNum exceeds the number of sequences in
// the file.
bool SequenceReaderMulti::findSequence( SequenceNumber seqNum )
{

  //  SequenceNumber currentSeqNum(seqNum);
  currentSeqNum_ = seqNum;
  for ( thisReader_ = allReaders_.begin(); 
        thisReader_ != allReaders_.end(); thisReader_++ )
  {
    if (thisReader_->allSeqsRead_)
    {
      if ( currentSeqNum_ <= thisReader_->size_ )
      {
        //      assert( thisReader_->ptr_->findSequence( currentSeqNum_ ) == true );
        lastSequenceNumber_ = --seqNum; // last read = 1 behind current
        return true;    
      } // ~if
    } // ~if
    else
    {
      if ( thisReader_->ptr_->findSequence( currentSeqNum_ ) == true )
      {
        lastSequenceNumber_ = --seqNum; // last read = 1 behind current
        return true;
      }
      else
        thisReader_->size_ = thisReader_->ptr_->getNumSequencesInFile();
      //      thisReader_->ptr_->rewind();
    } // ~else
    
    
    currentSeqNum_ -= thisReader_->ptr_->getNumSequencesInFile();


  } // ~for

  return false;

    //    if (thisReader_->ptr_->findSequence( currentSeqNum )) 
    //    {
      //      if (!thisReader_->allSeqsRead_)
      //   {
      //      thisReader_->allSeqsRead_=true;
      //      thisReader_->size_ = thisReader_->ptr_->getNumSequencesInFile();
      //    } // ~if
    //      return true;
    //   }
    //    if (!thisReader_->allSeqsRead_)
    //    {
    //     thisReader_->allSeqsRead_=true;
    //    thisReader_->size_ = thisReader_->ptr_->getNumSequencesInFile();
    //    } // ~if
    //  } // ~for
  //  return false;

} // ~SequenceReaderMulti::findSequence( SequenceNumber seqNum )




  // Function Name: findReader
  // Arguments: SequenceNumber& (in/out)
  // Given an input sequence number, adjusts seqNum and thisReader_
  // so that the seqNum'th sequence of *this is obtained by passing
  // the adjusted value of seqNum to *thisReader_
  bool SequenceReaderMulti::findReader( SequenceNumber& seqNum )
  {
    DEBUG_L3( "SequenceReaderMulti::findReader" );

    // NB calling getNumSequencesInFile also makes sure numSeqs_ has
    // been set up

    if ( seqNum > getNumSequencesInFile() )
    {
      return false;
      //      monitoringStream_ 
      //   << "findReader: Requested sequence number (" << seqNum 
      //   << ")\n is outside range of sequences in file (1 to " 
      //  << getNumSequencesInFile() << ")." << endl;
      //   throw SSAHAException("Invalid sequence number.");
    } // ~if

    for ( vector<SeqReaderInfo>::iterator i = allReaders_.begin(); 
          i != allReaders_.end(); i++ )
    {
      if ( seqNum <= i->size_ ) { thisReader_ = i; break; } // %%%%%
      seqNum -= i->size_;
    } // ~for i

    return true;
    
    // cout << allReaders_.begin().size_ << " " << seqNum << " " << i << endl;
  } // ~void SequenceReaderMulti::findReader( int& seqNum )

  // Function Name: changeMode
  // Arguments: const SequenceReaderMode&
  // Makes a copy of mode and uses it to handle mismatch character reads
  void SequenceReaderMulti::changeMode( SequenceReaderMode* pMode )
  {
    for ( vector<SeqReaderInfo>::iterator i = allReaders_.begin(); 
          i != allReaders_.end(); i++ ) 
          {
            (i->ptr_)->changeMode(pMode);
          }
  } // ~SequenceReaderMulti::changeMode( SequenceReaderMode* pMode )




  // Function Name: getNextSequence
  // Arguments: WordSequence& (out), int (in)
  // Returns:   int
  // Read the next set of sequence information from the file and parse it
  // into WordSequence format. Returns -1 if there has been a problem with
  // reading the sequence, else returns the number of valid base pairs 
  // contained within the final word of the sequence
  int SequenceReaderMulti::getNextSequence
  ( WordSequence& nextSeq, int wordLength )
  {
    DEBUG_L2( "SequenceReaderMulti::getNextSequence" );

    int numInLast;

    while 
      (    ( thisReader_ 
             != allReaders_.end() )
        && ( (   numInLast
               = thisReader_->ptr_->getNextSequence( nextSeq, wordLength ) )
              == -1 )    )
    {
      if ( thisReader_->allSeqsRead_ == false )
      {
        thisReader_->size_ = (thisReader_->ptr_)->getNumSequencesInFile();
        thisReader_->allSeqsRead_ = true;
      } // ~if
      thisReader_++;
      //      if (thisReader_!=allReaders_.end()) thisReader_->ptr_->rewind();
    } // ~while

    if (thisReader_ == allReaders_.end()) 
    {
      monitoringStream_ 
      << "SequenceReaderMulti::getNextSequence - last seq has been read.\n";
      return -1;
    } // ~if

    if ( numInLast != -1 ) lastSequenceNumber_++;

    return numInLast; 

  } // ~SequenceReaderMulti::getNextSequence


  // Function Name: getSequence
  // Arguments: WordSequence& (out), SequenceNumber (in), int (in)
  // Returns:   bool
  // Read the sequenceNumber-th set of sequence information from the file and 
  // parse it into WordSequence format
  int SequenceReaderMulti::getSequence
  ( WordSequence& nextSeq, SequenceNumber sequenceNumber, int wordLength )
  {
    DEBUG_L2( "SequenceReaderMulti::getSequence" );

    //    findReader( sequenceNumber ); 
    if (!findReader( sequenceNumber )) return -1;

    lastSequenceNumber_ = sequenceNumber;
            (thisReader_->ptr_)->rewind();

    return (thisReader_->ptr_)->getSequence
                                  ( nextSeq, sequenceNumber, wordLength );

  } // ~SequenceReaderMulti::getSequence

  // Function Name: getLastSequenceName
  // Arguments: string& (out)
  // Returns:   void
  // Fills the string with the name of the last sequence read 
  void SequenceReaderMulti::getLastSequenceName( string& seqName ) const 
  {
    (thisReader_->ptr_)->getLastSequenceName( seqName );
  } // ~SequenceReaderMulti::getLastSequenceName( string& seqName ) const 

  // Function Name: computeNumSequencesInFile
  // Arguments: void
  // Returns:   SequenceNumber
  // Returns the number of sequences in the file (will be done by lazy 
  // initialization, i.e. will only be calculated if asked for. NB this
  // will lose the current place in the file)
  SequenceNumber SequenceReaderMulti::computeNumSequencesInFile( void )
  {
    SequenceNumber numSeqs = 0;
    for ( vector<SeqReaderInfo>::iterator i 
          = allReaders_.begin(); 
          i != allReaders_.end();
          i++ )
    {
      if ( i->allSeqsRead_ == false )
      {
        i->size_ =  (i->ptr_)->getNumSequencesInFile();
        i->allSeqsRead_ = true;
      } // ~if
      numSeqs += i->size_;       
    } // ~for i
    thisReader_ = allReaders_.end();
    lastSequenceNumber_ = numSeqs;
    return numSeqs; 
  } // ~SequenceReaderMulti::computeNumSequencesInFile( void )

  // Function Name: getBitsPerSymbol
  // Arguments: none
  // Returns:   int
  // Returns number of bits per symbol used in encoding
int SequenceReaderMulti::getBitsPerSymbol ( void ) const
{
  if (isFirstSeq_) assert(1==0);
  return bitsPerSymbol_;
  // fails if thisReader is at end 
  //  return ((thisReader_->ptr_)->getBitsPerSymbol());
} // ~SequenceReaderMulti::getBitsPerSymbol ( void ) const

// Function Name: getSourceDataType
// Arguments: none
// Returns:   SourceDataType
// Returns type of data being encoded (protein or DNA)
SourceDataType SequenceReaderMulti::getSourceDataType( void ) const
{
  if (isFirstSeq_) assert(1==0);
  return sourceDataType_;
  // fails if thisReader is at end 
  //  return ((thisReader_->ptr_)->getSourceDataType());
} // ~SequenceReaderMulti::getSourceDataType ( void ) const

  // Function Name: printName
  // Arguments: string& (out), SequenceNumber (in)
  // Returns:   void
  // Fills a string with the name of the requested sequence
bool SequenceReaderMulti::printName( ostream& os, SequenceNumber seqNum ) 
{

  lastSequenceNumber_ = seqNum;
  //  findReader( seqNum ); 
  if (!findReader(seqNum)) return false;
  (thisReader_->ptr_)->printName( os, seqNum );
  return true;
} // ~SequenceReaderMulti::printName

  // Function Name: printSideInfo
  // Arguments: string& (out), SequenceNumber (in)
  // Returns:   void
  // Fills a string with the name of the requested sequence
bool SequenceReaderMulti::printSideInfo( ostream& os, SequenceNumber seqNum )
{

  lastSequenceNumber_ = seqNum;
  //  findReader( seqNum ); 
  if (!findReader(seqNum)) return false;
  (thisReader_->ptr_)->printSideInfo( os, seqNum );
  return true;
} // ~SequenceReaderMulti::printName


  // Function Name: printSource
  // Arguments: string& (out), SequenceNumber (in)
  // Returns:   void
  // Fills a string with the name of the requested sequence
bool SequenceReaderMulti::printSource( ostream& os, SequenceNumber seqNum )
{
  lastSequenceNumber_ = seqNum;
  //findReader( seqNum ); 
  if (!findReader(seqNum)) return false;
  (thisReader_->ptr_)->printSource( os, seqNum );
  return true;
} // ~SequenceReaderMulti::printName

// Function Name: extractSource
// This extracts the source data for bases seqStart to seqEnd inclusive
// of sequence seqNum and places it in source
void SequenceReaderMulti::extractSource
( char** pSource,//vector<char>& source, 
  SequenceNumber seqNum,
  SequenceOffset seqStart,
  SequenceOffset seqEnd )
{
  //  cout << "SRM::extractSource " << seqNum << endl;
  //  SequenceNumber savedSeqNum(lastSequenceNumber_);
  //  vector<SeqReaderInfo>::iterator savedReader(thisReader_);
  //  SequenceNumber savedSeqNumPtr(thisReader_->ptr_->getLastSequenceNumber());
  //  monitoringStream_ << "srm::Extract1 " << seqNum << " " << thisReader_ << endl;

  SequenceReaderState* pState(saveState());

  if (!findSequence(seqNum))
  { 
    monitoringStream_ << "extSource: Requested sequence number (" << seqNum
                      << ") exceeds number of sequences in collection ("
                      << getNumSequencesInFile() << ")." << endl;
    throw SSAHAException
      ("Invalid sequence number in SequenceReaderMulti::extractSource");
  } // ~if   
  //  monitoringStream_ << "srm::Extract2 " << seqNum << " " << thisReader_ << endl;
  (thisReader_->ptr_)->extractSource
    ( pSource, currentSeqNum_, seqStart, seqEnd );

  //  lastSequenceNumber_ = savedSeqNum;

  //  thisReader_->ptr_->rewind();
  //  findSequence( savedSeqNum+1 );

  restoreState(pState);
  //  lastSequenceNumber_=savedSeqNum;
  //  thisReader_=savedReader;
  //  thisReader_->ptr_->findSequence(savedSeqNumPtr+1);

} // ~SequenceReaderMulti::extractSource


void SequenceReaderMulti::saveIndexImp
( ostream& fileFile, 
  ostream& indexFile, 
  int& fileNumber )
{
  computeNumSequencesInFile();
  for ( vector<SeqReaderInfo>::iterator i 
          = allReaders_.begin(); 
          i != allReaders_.end();
          i++ )
  {
    (i->ptr_)->saveIndexImp( fileFile, indexFile, fileNumber );
    fileNumber++;
  } // ~for i

} // ~SequenceReaderMulti::saveIndexImp



// End of file SequenceReaderMulti.cpp

---