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

// Description:

// Includes:

#include "SequenceReader.h"
#include "HashTable.h"
#include <string>
#include <fstream>

static const char nc='\0';
static const char reverseChar[256] = 
{
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,'-',nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc, /* next is 'A' */ 
'T',nc,'G',nc,nc,nc,'C',nc,nc,nc,nc,nc,nc,'N',nc, /* next is 'P' */
nc,nc,nc,nc,'A',nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc, /* next is 'a' */ 
't',nc,'g',nc,nc,nc,'c',nc,nc,nc,nc,nc,nc,'n',nc, /* next is 'p' */
nc,nc,nc,nc,'a',nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,
nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc,nc
};





// ### Function Definitions ###

SequenceReader::SequenceReader
( ostream& monitoringStream  ) : 
  //  pInputFileStream_( new ifstream( fileName.c_str() ) ), 
  //  inputFileStream_( fileName.c_str() ), 
getName( this  ),
getSideInfo( this ),
getSource( this ),
lastSequenceNumber_(0),
allSequencesRead_(false),
numSequencesInFile_(0),
monitoringStream_( monitoringStream )
{
  monitoringStream_ << "constructing SequenceReader" << endl;
} // constructor

SequenceReader::SequenceReader( const SequenceReader& rhs ) :
  getName( rhs.getName ),
  getSideInfo( rhs.getSideInfo ),
  getSource( rhs.getSource ),
  lastSequenceNumber_( rhs.lastSequenceNumber_ ),
  allSequencesRead_( rhs.allSequencesRead_ ),
  numSequencesInFile_( rhs.numSequencesInFile_ ),
  monitoringStream_(   rhs.monitoringStream_ )
{
  monitoringStream_ << "copy constructing SequenceReader" << endl;
} // copy constructor

SequenceReader::~SequenceReader()
{
  monitoringStream_ << "destructing SequenceReader" << endl;
} // destructor

  // Function Name: changeMode
  // Arguments: const SequenceReaderMode&
  // Makes a copy of mode and uses it to handle mismatch character reads
// Now pure virtual in this class TC 14.3.01
//  void SequenceReader::changeMode( SequenceReaderMode* pMode )
//  {
//    monitoringStream_ << "SequenceReader::changeMode\n";
//    delete pState_;
//    pState_ = pMode->clone();
//  }


  // Function Name: getSequence
  // Arguments: WordSequence& (out), SequenceNumber (in), const HashTable& (in)
  // Returns:   bool
  // Read the sequenceNumber-th set of sequence information from the file and 
  // parse it into WordSequence format, getting word length from hashTable
  int SequenceReader::getSequence
  ( WordSequence& nextSeq, SequenceNumber sequenceNumber, const HashTable& hashTable )
  {
      return getSequence( nextSeq, sequenceNumber, hashTable.getWordLength() );
  } // ~getSequence

  // Function Name: getNextSequence
  // Arguments: WordSequence& (out), HashTable& (in)
  // Returns:   int
  // Read the next set of sequence information from the file and parse it
  // into WordSequence format, using word length required by hashTable
  int SequenceReader::getNextSequence( WordSequence& nextSeq, const HashTable& hashTable )
  {
      return getNextSequence( nextSeq, hashTable.getWordLength() );
  } // ~getNextSequence


void SequenceReaderNamePrinter::print( ostream& os )
{ 
  if (!pReader_->printName(os,seqNum_)) throw NumberOutOfRange();
} // ~print

void SequenceReaderSideInfoPrinter::print( ostream& os )
{ 
  if (!pReader_->printSideInfo(os,seqNum_)) throw NumberOutOfRange();
} // ~print

void SequenceReaderSourcePrinter::print( ostream& os )
{ 
  if (!pReader_->printSource(os,seqNum_)) throw NumberOutOfRange();
} // ~print

// SourceReader member function definitions

void SourceReader::extractSource
( char** pSource, //vector<char>& source, 
  SequenceNumber seqNum,
  SequenceOffset seqStart,
  SequenceOffset seqEnd )
{
  throw SSAHAException
    ("SourceReader::extractSource not implemented for this subclass!");
} // ~SourceReader::extractSource

void SourceReader::extractSourceReverse
( char** pSource, //vector<char>& source, 
  SequenceNumber seqNum,
  SequenceOffset seqStart,
  SequenceOffset seqEnd )
{
  //  vector<char> temp;
  //  extractSource( temp, seqNum, seqStart, seqEnd );
  //  source.resize( source.size() + temp.size() );
  //  vector<char>::iterator i(temp.begin()), j(&source.back());

  reverseBuffer_.resize( seqEnd - seqStart + 1 );

  char* pSourceFwd;
  extractSource( &pSourceFwd, seqNum, seqStart, seqEnd );

  
  vector<char>::iterator i(reverseBuffer_.begin());
  char* j = &pSourceFwd[ seqEnd - seqStart ];

  for ( ; i != reverseBuffer_.end() ; i++, j-- )
  {
    *i= reverseChar[ *j ];
    //    cout << "char: " << (int)*j << *j << " - " << (int)*i << *i << endl;
  }

  *pSource = (char*)&reverseBuffer_[0];

} // ~SourceReader::extractSourceReverse







// save the location of the start of each sequence in a file
// sets up ofstream and calls the virtual saveIndex below
void SourceReader::saveIndex( const string& fileName )
{
  ofstream filesFile((fileName+(string)".files").c_str());
  ofstream indexFile((fileName+(string)".index").c_str());
  assert(!filesFile.fail());
  assert(!indexFile.fail());
  int fileNum(0);
  saveIndexImp(filesFile, indexFile, fileNum);
  filesFile.close();
  indexFile.close();
} // ~SourceReader::saveIndex 


void SourceReader::saveIndexImp
( ostream& fileFile, 
  ostream& indexFile, 
  int& fileNumber )
{
  throw SSAHAException("SourceReader::saveIndex not implemented for this subclass!");
}

SourceReaderIndex::SourceReaderIndex( const string& fileName ) :
currentFileNum_(0)
{
  ifstream filesFile((fileName+(string)".files").c_str());
  ifstream indexFile((fileName+(string)".index").c_str());
  assert(!filesFile.fail());
  assert(!indexFile.fail());

  // read in file names
  fileNames_.push_back( new string );
  
  while( filesFile >> *(fileNames_.back()) ) 
    fileNames_.push_back( new string );

  // get rid of unwanted empty string at end
  delete fileNames_.back();
  fileNames_.pop_back();

  filesFile.close();
  // read in index

  indexFile.seekg(0, ios::end);
  std::streampos indexSize( indexFile.tellg() );
  indexFile.seekg(0, ios::beg);

  assert ( indexSize % sizeof (SeqIndexInfo) == 0 );

  numSeqs_ = indexSize / sizeof( SeqIndexInfo );

  index_.resize( numSeqs_ );

  indexFile.read( (char*)&index_[0], indexSize );

  indexFile.close();
 
  //  for ( int i(0) ; i < index_.size() ; i++ )
  //    cout << "index: " << i << " " << index_[i].fileNum << " " 
  // << index_[i].seqPos << endl;


  pCurrentFile_ =  new ifstream( (*fileNames_[0]).c_str() );


} // ~SourceReaderIndex::SourceReaderIndex( const string& fileName )


SourceReaderIndex::~SourceReaderIndex()
{
  for ( vector<string*>::iterator i( fileNames_.begin() );
       i != fileNames_.end() ; ++i ) delete *i;
} // ~SourceReaderIndex()::~SourceReaderIndex()

const char* SourceReaderIndex::extractName( SequenceNumber seqNum )
{
  if (seqNum != lastNameSeqNum_)
  { // need to read a new sequence name

    assert( seqNum !=0 );
    assert( seqNum <= index_.size() );

    vector<SeqIndexInfo>::iterator pThisSeq(&index_[seqNum-1]);

    // cout << "seq location " << pThisSeq->fileNum << " " << pThisSeq->seqPos << endl;
    // cout << *fileNames_[pThisSeq->fileNum] << endl;

    if (pThisSeq->fileNum != currentFileNum_ )
    { // need to move to a different sequence file
      pCurrentFile_->close();
      delete pCurrentFile_;
      pCurrentFile_ 
        =  new ifstream( (*fileNames_[pThisSeq->fileNum]).c_str() );
      assert(!pCurrentFile_->fail());
      currentFileNum_ = pThisSeq->fileNum;
    }

    // wind to start of sequence
    pCurrentFile_->clear();
    pCurrentFile_->seekg( pThisSeq->seqPos, ios::beg );

    pCurrentFile_->getline
      ( nameBuffer_, nameBufferSize_, '\n' );
    lastNameSeqNum_ = seqNum; // +1 cos we decremented seqNum earlier
    char* pSpace = strchr( nameBuffer_, ' ' );
    if (pSpace!=NULL) *pSpace='\0';
    
  } // ~if
  assert((nameBuffer_[0]=='>')||(nameBuffer_[0]=='@'));
  return &nameBuffer_[1];

} // ~SourceReaderIndex::extractName( SequenceNumber seqNum )



void SourceReaderIndex::extractSource
( char** pSource, //vector<char>& source, 
  SequenceNumber seqNum,
  SequenceOffset seqStart,
  SequenceOffset seqEnd )
{

  if (seqNum != lastSourceSeqNum_)
  { // need to put a new sequence in the cache

    assert( seqNum !=0 );
    assert( seqNum <= index_.size() );

    vector<SeqIndexInfo>::iterator pThisSeq(&index_[seqNum-1]);

    // cout << "seq location " << pThisSeq->fileNum << " " << pThisSeq->seqPos << endl;
    // cout << *fileNames_[pThisSeq->fileNum] << endl;

    if (pThisSeq->fileNum != currentFileNum_ )
    { // need to move to a different sequence file
      pCurrentFile_->close();
      delete pCurrentFile_;
      pCurrentFile_ 
        =  new ifstream( (*fileNames_[pThisSeq->fileNum]).c_str() );
      assert(!pCurrentFile_->fail());
      currentFileNum_ = pThisSeq->fileNum;
    }

    // wind to start of sequence
    pCurrentFile_->clear();
    pCurrentFile_->seekg( pThisSeq->seqPos, ios::beg );
    //  cout << "about to extract to cache "<< pCurrentFile_->tellg() << endl;
    extractToCache( pCurrentFile_ );
    lastSourceSeqNum_ = seqNum; // +1 cos we decremented seqNum earlier
    //    cout << "cached seq " << lastSourceSeqNum_ << " " << lastSourceSeq_.size() << endl;
    lastNameSeqNum_ = seqNum;
    char* pSpace = strchr( nameBuffer_, ' ' );
    if (pSpace!=NULL) *pSpace='\0';
  } // ~if

  // check for japes and tomfoolery

  if ( seqStart>seqEnd)
  {
      throw SSAHAException
        ("Requested seq start exceeds requested seq end in SourceReaderIndex::extractSource");
  } // ~if
  else if (seqEnd>lastSourceSeq_.size() )
  {
    cout << seqEnd << " " << lastSourceSeq_.size() << endl;
      throw SSAHAException
        ("Requested last byte exceeds end of seq in SourceReaderIndex::extractSource");
  } // ~else if


  *pSource = &lastSourceSeq_[seqStart-1]; 

  //  vector<char>::size_type currentSize(source.size());
  //  source.resize(currentSize+seqEnd-seqStart+1);

  //  memcpy( (char*)&source[currentSize], (char*)&lastSourceSeq_[seqStart-1], seqEnd-seqStart+1);


}

void SourceReader::extractToCache( istream* pCurrentFile)
{
   // Check that this is a sensible place for sequence start
    int firstOfLine(pCurrentFile->peek());
    assert(((char)firstOfLine=='>')||((char)firstOfLine=='@'));

    int currentSize;


    lastSourceSeq_.resize(sourceBufferSize_);

    char* pStart = (char*) &lastSourceSeq_[0];
    char* pResize = (char*) &lastSourceSeq_.back()-resizeCacheThreshold_;

    //  get sequence name and put in nameBuffer_
    //    pCurrentFile->getline
    //  ( (char*)&lastSourceSeq_[0], sourceBufferSize_, '\n' );
    pCurrentFile->getline
      ( nameBuffer_, nameBufferSize_, '\n' );

    std::streampos lastPos;
    //    int firstOfLine;

    while(1)
    {
      lastPos = pCurrentFile->tellg();
      firstOfLine = pCurrentFile->peek();
      if ( ((char)firstOfLine=='>') || 
           ((char)firstOfLine=='+') || 
           (firstOfLine==EOF) ) break;
      pCurrentFile->getline( pStart, sourceBufferSize_, '\n' );
      pStart += (unsigned long)pCurrentFile->tellg()-lastPos-1; // -1 accounts for \n at end
      if (pStart >= pResize)
      {
        currentSize = pStart-&lastSourceSeq_[0]; 
        // need currentSize because storage for lastSourceSeq_ may move
        // when resized, thus invalidating pStart
        lastSourceSeq_.resize
          ( lastSourceSeq_.size() + (lastSourceSeq_.size()>>1) );
        pStart = &lastSourceSeq_[currentSize];
        pResize = (char*) &lastSourceSeq_.back()-resizeCacheThreshold_;
      } // ~if
    } // ~while

    lastSourceSeq_.resize( pStart-&lastSourceSeq_[0] );
    //    cout << "Read " << lastSourceSeq_.size() << " characters into cache.\n";



}






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

// End of file SequenceReader.cpp

---