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

// Description:

// Includes:

#include "SequenceEncoder.h"
#include "SequenceReader.h"
#include "SequenceReaderFasta.h"
#include <strstream>
#include <fstream>


// ### Function Definitions ###

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

  SequenceReaderFasta::SequenceReaderFasta
  ( const char* fileName, SequenceEncoder* pEncoder,
    ostream& monitoringStream )
   : SequenceReaderFile ( fileName, '>', '>', pEncoder->clone(),
                         monitoringStream )
  {
    monitoringStream_ << "constructing SequenceReaderFasta"  << endl;
  } // ~constructor

  SequenceReaderFasta::SequenceReaderFasta
  ( const char* fileName,
    ostream& monitoringStream )
   : SequenceReaderFile ( fileName, '>', '>', new SequenceEncoderDNA(12),
                         monitoringStream )
  {
    monitoringStream_ << "constructing SequenceReaderFasta" << endl;
  } // ~constructor

  SequenceReaderFastaProtein::SequenceReaderFastaProtein
  ( const char* fileName, ostream& monitoringStream ) 
  : SequenceReaderFasta ( fileName, new SequenceEncoderProtein(5), 
                          monitoringStream )
  {
    monitoringStream_ << "constructing SequenceReaderFastaProtein" << endl;
  } // ~constructor

  SequenceReaderFile::SequenceReaderFile
  ( const char* fileName,
    char seqStartChar,
    char seqStopChar,
    SequenceEncoder* pEncoder,
    ostream& monitoringStream )
    : SequenceReader( monitoringStream ),
      seqStartChar_( seqStartChar ),
      seqStopChar_( seqStopChar ),
      //      lastSourceSeqNum_( 0 ),
      pInputFileStream_( new ifstream( fileName ) ),
      fileName_( fileName ),
      pEncoder_( pEncoder )
  {

    monitoringStream_ << "constructing SequenceReaderFile" << this << endl;
    if ( pInputFileStream_->fail() )
    {
      throw SSAHAException
      (   (string)"SequenceReaderFile - unable to open input file "
        + (string)fileName );

    } // ~if
    monitoringStream_ << "SequenceReaderFile::connected to file "
                      << fileName << endl;

  } // ~constructor

  SequenceReaderFile::SequenceReaderFile
  ( istream& inputStream,
    char seqStartChar,
    char seqStopChar,
    SequenceEncoder* pEncoder,
    ostream& monitoringStream )
    : SequenceReader( monitoringStream ),
      seqStartChar_( seqStartChar ),
      seqStopChar_( seqStopChar ),
      //   lastSourceSeqNum_( 0 ),
      pInputFileStream_( &inputStream ),
      //      fileName_( fileName ),
      pEncoder_( pEncoder )
  {

    monitoringStream_ << "constructing SequenceReaderFile" << this << endl;
    if ( pInputFileStream_->fail() )
    {
      throw SSAHAException
      ( "SequenceReaderFile - unable to open input stream " );

    } // ~if
    monitoringStream_ << "SequenceReaderFile::connected to input file stream."
                      << endl;

  } // ~constructor



    

  SequenceReaderFile::SequenceReaderFile( const SequenceReaderFile& rhs) 
  : SequenceReader( rhs.monitoringStream_ ),
    seqStartChar_( rhs.seqStartChar_ ),
    seqStopChar_( rhs.seqStopChar_ ),
    pInputFileStream_( new ifstream( rhs.fileName_.c_str() ) ),
    fileName_( rhs.fileName_.c_str() ),
    seqPositions_( rhs.seqPositions_ ),
    //   lastSourceSeqNum_(0),
    pEncoder_( rhs.pEncoder_->clone() )
  {
    monitoringStream_ << "copy constructing SequenceReaderFile" << this 
                      << endl;
    if (pInputFileStream_->fail())
    {
       throw SSAHAException
             (   (string)"SequenceReaderFile - unable to open input file " 
               + (string)fileName_ );

    } // ~if

    monitoringStream_ << "SequenceReaderFile:SequenceReader: "
                      << "connected to file " << fileName_ << endl;

  } // ~copy constructor

  SequenceReaderFile::~SequenceReaderFile()
  {
    monitoringStream_ << "destructing SequenceReaderFile" << this << endl;
    delete pEncoder_;
    ifstream*  pf(dynamic_cast<ifstream*>(pInputFileStream_));
    if (pf!=NULL) 
    {
        pf->close();
        delete pInputFileStream_;
    }
  } // ~destructor

  // Function Name: changeMode
  // Arguments: const SequenceReaderMode&
  // Makes a copy of mode and uses it to handle mismatch character reads
  void SequenceReaderFile::changeMode( SequenceReaderMode* pMode )
  {
    pEncoder_->changeMode( pMode );
  }

  // Function Name: rewind
  // Arguments: void
  // Returns:   void
  // Rewind to the start of the data file, so that getNextSequence will
  // return the first sequence in the file. NB Also clears any error flags
  // that are set on pInputFileStream_ 
  void SequenceReaderFile::rewind( void ) 
  {
    // rewind to start of file
    pInputFileStream_->clear();
    pInputFileStream_->seekg( 0 );
    lastSequenceNumber_ = 0;
  } // ~SequenceReaderFile::rewind( void )

  // Function Name: findSequence
  // Arguments: SequenceNumber (in)
  // Returns:   void
  // Winds the input file stream to the start of sequence number seqNum. 
  // Throws an exception if seqNum exceeds the number of sequences in
  // the file.
  bool SequenceReaderFile::findSequence( SequenceNumber seqNum )
  {
    
    //    cout << "SRF::findSeq start" << pInputFileStream_->tellg() << endl; 
    pInputFileStream_->clear();
    if ( seqNum ==0 )
    { 
      monitoringStream_ 
        << "SRF::findSeq: Requested sequence number (0) is not valid (must start at 1)." 
        << endl;
      return false;
    } // ~if
    else if ( seqNum <= seqPositions_.size() )
    {
      pInputFileStream_->seekg(seqPositions_[ seqNum - 1 ],ios::beg);
    } // ~if
    else if ( allSequencesRead_ == false )
    { // need to find more sequences in the file

      if ( seqPositions_.size() != 0 )
      {
        pInputFileStream_->seekg
        (seqPositions_.back(), ios::beg);
        pInputFileStream_->getline
        ( inputBuffer_ , inputBufferSize_, '\n' );
      } // ~if
   
      while(1)
      {
        //        pInputFileStream_->getline
        //      ( inputBuffer_ , inputBufferSize_, '\n' );
              //     getline( *pInputFileStream_, inputBuffer_, '\n' );
        if (pInputFileStream_->peek() == EOF) 
        {
          numSequencesInFile_ = seqPositions_.size();
          //      cout << "SRF::findSeq: asr = true" << endl;
          allSequencesRead_ = true;
          monitoringStream_ 
          << "SRF::findSeq(EOF): Requested sequence number (" << seqNum 
          << ")\nis outside range of sequences in file (1 to " 
          << numSequencesInFile_ << ")." << endl;
          //          throw NumberOutOfRange();
          // cout << "SRF::findSeq end" << pInputFileStream_->tellg() << endl; 
          return false;
        } // ~if
        else if (pInputFileStream_->peek() == seqStartChar_) 
        {
          seqPositions_.push_back(pInputFileStream_->tellg());
          if( seqPositions_.size() == seqNum ) break;
        } // ~else if
        pInputFileStream_->getline
        ( inputBuffer_ , inputBufferSize_, '\n' );
      } // ~while
      //      while ( seqPositions_.size() != seqNum );

    lastSequenceNumber_ = seqNum - 1;
    
    } // ~else if
    else
    { 
      monitoringStream_ 
        << "SRF::findSeq(asr): Requested sequence number (" << seqNum 
        << ")\nis outside range of sequences in file (1 to " 
        << numSequencesInFile_ << ")." << endl;
      //      throw NumberOutOfRange();

      //      if (seqNum == numSequencesInFile_+1)
      //     {
      //        cout << "SRF::findSeq - flicking to end" << endl; 
      //        pInputFileStream_->seekg( 0, ios::end );
      //  }
      //  cout << "SRF::findSeq end" << pInputFileStream_->tellg() << endl; 
      return false;
    } // ~else

    lastSequenceNumber_ = seqNum - 1;
    //  cout << "SRF::findSeq end" << pInputFileStream_->tellg() << endl; 
    return true;
  } // ~SequenceReaderFile::findSequence( SequenceNumber seqNum )


  // Function Name: getNextSequence
  // Arguments: WordSequence& (out), int (in)
  // Returns:   bool
  // Read the next set of sequence information from the file and parse it
  // into WordSequence format
  int SequenceReaderFile::getNextSequence
  ( WordSequence& nextSeq, int wordLength )
  {
    DEBUG_L2( "SequenceReaderFile::getNextSequence" );
    //   cout << "SRF::getNext start" << pInputFileStream_->tellg() << endl; 
    //    char firstOfLine;
    // `Interesting' standard library quirk: even though we are exclusively
    // reading chars firstOfLine must be an int (cos EOF is an int not a char)
    int firstOfLine;
    pInputFileStream_->clear();
    firstOfLine = pInputFileStream_->peek();
    if ( firstOfLine == EOF )
    {
      monitoringStream_ << "End of file has been reached." << endl;
      if ( allSequencesRead_ == false )
      {
        monitoringStream_ << "Setting numSequencesInFile to " 
                          << lastSequenceNumber_ << "." << endl;
        numSequencesInFile_ = lastSequenceNumber_;
        //  cout << "SRF::getNextSeq: asr = true" << endl;
        //  cout << pInputFileStream_->tellg() << endl;
          assert(!pInputFileStream_->fail());
          
        allSequencesRead_ = true;
      } // ~if
      // cout << "SRF::getNext end" << pInputFileStream_->tellg() << endl; 
      return -1;
    } // ~if

    std::streampos startPos( pInputFileStream_->tellg() );     

    if ( (char)firstOfLine == seqStartChar_ )
    { // then everything until the next '\n' is side info 
      //      pInputFileStream_->getline( sideInfoBuffer_, sideInfoBufferSize_, '\n' );
      getline( *pInputFileStream_, sideInfoBuffer_, '\n' );
      DEBUG_L2( "Read " << pInputFileStream_->gcount()
                        << " bytes of side info to buffer" );
    } // ~if
    else
    {
      monitoringStream_ 
        << "Error: file not in expected format (expected \"" 
        << seqStartChar_ << "\", received \"" << (char) firstOfLine
        << "\" instead).\n";
      throw SSAHAException("Data in wrong format!");
      //      sideInfoBuffer_="No side info for sequence.";
    } // ~else

    //    int basesInSequence(0); // this accumulates to the total number of
    // bases in the sequence (may be spread across multiple lines)
    //    Word thisWord(0);

    pEncoder_->setWordLength(wordLength);
    pEncoder_->linkSeq(nextSeq);

    while (1==1)
    {

      int numChars; // this is the number of bases on a line of FASTA file
                    // info (i.e. from one '\n' to the next) 
      firstOfLine = pInputFileStream_->peek();

      if (((char)firstOfLine == seqStopChar_) || ( firstOfLine == EOF )) break;
      pInputFileStream_->getline( inputBuffer_, inputBufferSize_, '\n' );
      //          getline( *pInputFileStream_, inputBuffer_, '\n' );

      DEBUG_L2( "Read " << pInputFileStream_->gcount()
                        << " bytes of sequence data to buffer" );


          numChars = pInputFileStream_->gcount() - 1; 
          // the '-1' accounts for '\n' at end of read
          //          numChars = inputBuffer_.size();

      pEncoder_->encode( inputBuffer_, numChars );

    }; // ~while

    pEncoder_->unlinkSeq();

    //    shiftLastWord( nextSeq, basesInSequence, wordLength );

    // If the last word is not completely full (as will be the case if
    // wordLength does not divide numBases) then shift the word so the
    // valid bases occupy the most significant bits of the word.
    //    nextSeq.numBasesInLast = basesInSequence % wordLength;
    //    thisWord = ( numBasesInLast == 0 ) ? 0 
    //     : thisWord << 2 * ( wordLength - numBasesInLast - 1);
    //   nextSeq.push_back(thisWord);
    //    if ( nextSeq.numBasesInLast != 0 )  
    //    nextSeq.back() <<= 2 * ( wordLength - nextSeq.numBasesInLast - 1);

    //    DEBUG_L2(    "There were " << basesInSequence
    //        << " base pairs in this sequence." );
      
    //    DEBUG_L2(    "Last word (" << printWord( nextSeq.back(), wordLength )
    //        << ") contains " << nextSeq.numBasesInLast
    //        << " valid base pairs" );

    DEBUG_L2("Finished reading sequence " << lastSequenceNumber_ );

    lastSequenceNumber_++;

    if ( lastSequenceNumber_ > seqPositions_.size() )
    {
      seqPositions_.push_back( startPos );
    } // ~if
    //    cout << "SRF::getNext end" << pInputFileStream_->tellg() << endl; 

    return nextSeq.getNumBasesInLast(); 

  } // ~SequenceReaderFile::getNextSequence

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

    if (!findSequence( sequenceNumber )) return -1;
    else return getNextSequence( nextSeq, wordLength );

  } // ~SequenceReaderFile::getSequence 

  // Function Name: getLastSequenceName
  // Arguments: string& (out)
  // Returns:   void
  // Fills the string with the name of the last sequence read 
  void SequenceReaderFile::getLastSequenceName( string& seqName ) const
  {
    DEBUG_L2( "SequenceReaderFile::getLastSequenceName" );
    string::size_type nameEnd = sideInfoBuffer_.find_first_of(' ');
    if ( nameEnd == string::npos ) nameEnd = sideInfoBuffer_.size();
    seqName = sideInfoBuffer_.substr( 1, nameEnd-1 );
    DEBUG_L2( "Last sequence name: " << seqName );
  } // ~SequenceReaderFile::getLastSequenceName

  // Function Name: printName
  // Arguments: ostream& (out), SequenceNumber (in)
  // Returns:   void
  // Sends the name of the requested sequence to the output stream.
  bool SequenceReaderFile::printName
  ( ostream& os, SequenceNumber seqNum )
  {
    DEBUG_L3( "SequenceReaderFile::printName" );

    if (!findSequence( seqNum )) return false;


    string firstLine;

    getline(*pInputFileStream_,firstLine,'\n');

    string::size_type nameEnd = firstLine.find_first_of(' ');
    if ( nameEnd == string::npos ) nameEnd = firstLine.size();

    os << firstLine.substr( 1, nameEnd ); 

    findSequence( seqNum );

    return true;
  } // ~SequenceReaderFile::printName


  // Function Name: printSideInfo
  // Arguments: ostream& (out), SequenceNumber (in)
  // Returns:   void
  // Sends the side info (e.g. clone name) for the requested sequence 
  // to the output stream.
  bool SequenceReaderFile::printSideInfo
  ( ostream& os, SequenceNumber seqNum )
  {
    DEBUG_L2( "SequenceReaderFile::printSideInfo" );
 
    if (!findSequence( seqNum )) return false;

    string firstLine;

    getline(*pInputFileStream_,firstLine,'\n');

    string::size_type nameEnd = firstLine.find_first_of(' ');

    if ( nameEnd == string::npos ) os << "No side info for this sequence.";
    else os << firstLine.substr(firstLine.find_first_of(' ') ); 

    findSequence( seqNum );

    return true;
  } // ~SequenceReaderFile::printSideInfo

  // Function Name: printSource
  // Arguments: ostream& (out), SequenceNumber (in)
  // Returns:   void
  // Send to the output stream the source data (in general, ASCII) from which 
  // the requested sequence was decoded.
  bool SequenceReaderFile::printSource
  ( ostream& os, SequenceNumber seqNum )
  {
    DEBUG_L2( "SequenceReaderFile::getSequenceSource" );

    if (!findSequence( seqNum )) return false;

    int firstOfLine((int)'x');
    while( ( (char)firstOfLine != seqStartChar_ ) && ( firstOfLine != EOF ) )
    {
      pInputFileStream_->getline( inputBuffer_, inputBufferSize_, '\n' );
      os << inputBuffer_ << endl;
      firstOfLine = pInputFileStream_->peek();
    } // ~while

    findSequence( seqNum );

    return true;
  } // ~SequenceReaderFile::printSource

  // Function Name: computeNumSequencesInFile
  // Arguments: void
  // Returns:   int
  // 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 SequenceReaderFile::computeNumSequencesInFile( void ) 
  {
    DEBUG_L2( "SequenceReaderFile::computeNumSequencesInFile" );

    SequenceNumber numSeqs = 0;

    // rewind to start of file
    pInputFileStream_->clear();
    pInputFileStream_->seekg( 0 );
    seqPositions_.clear();

    while ( pInputFileStream_->peek() != EOF )
    {
      if (pInputFileStream_->peek() == seqStartChar_) 
      {
        numSeqs++;
        seqPositions_.push_back(pInputFileStream_->tellg());
      } // ~if
              pInputFileStream_->getline
               ( inputBuffer_ , inputBufferSize_, '\n' );
              //     getline( *pInputFileStream_, inputBuffer_, '\n' );
    } // ~while

    lastSequenceNumber_ = numSeqs;

    DEBUG_L2( "Found " << numSeqs << " sequences in the file." );

    return numSeqs;

  } // ~SequenceReaderFile::computeNumSequencesInFile( void ) const

  // Function Name: getBitsPerSymbol
  // Arguments: none
  // Returns:   int
  // Returns number of bits per symbol used in encoding
int SequenceReaderFile::getBitsPerSymbol ( void ) const
{
  return pEncoder_->getBitsPerSymbol();
} // ~SequenceReaderFile::getBitsPerSymbol ( void ) const

  // Function Name: getSourceDataType
  // Arguments: none
  // Returns:   SourceDataType
  // Returns type of data being encoded (protein or DNA)
SourceDataType SequenceReaderFile::getSourceDataType( void ) const
{
  return pEncoder_->getSourceDataType();
}


// Function Name: extractSource
// This extracts the source data for bases seqStart to seqEnd inclusive
// of sequence seqNum and places it in source
void SequenceReaderFile::extractSource
( char** pSource, //vector<char>& source, 
  SequenceNumber seqNum,
  SequenceOffset seqStart,
  SequenceOffset seqEnd )
{

  if (seqNum!= lastSourceSeqNum_)
  {

    SequenceReaderState* pState(saveState());

    if (!findSequence( seqNum )) 
    {

      throw SSAHAException
        ("Could not find start of seq in SequenceReaderFile::extractSource");
    } // ~if

    //    cout << "moved to " << pInputFileStream_->tellg() << endl;

    // Check that this is a sensible place for sequence start
    //    int firstOfLine(pInputFileStream_->peek());
    //  if (((char)firstOfLine)!=seqStartChar_)
    //   {
    //   throw SSAHAException
    //  ("Did not find expected start char in SequenceReaderFile::extractSource");
    //  } // ~if    


    extractToCache( pInputFileStream_ );
    lastSourceSeqNum_=seqNum;

    // return to original position in file
    restoreState(pState);

  } // ~if

  // check for japes and tomfoolery

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

  *pSource = &lastSourceSeq_[seqStart-1];
  // copy relevant bit of sequence across (TBD this is a waste, just return ptr)
  //  vector<char>::size_type currentSize(source.size());
  //  source.resize(currentSize+seqEnd-seqStart+1);
  //  memcpy( (char*)&source[currentSize], (char*)&lastSourceSeq_[seqStart-1], seqEnd-seqStart+1);
  //  cout << "SRF::extract end" << pInputFileStream_->tellg() << endl; 

} // ~SequenceReaderFile::extractSource



void SequenceReaderFile::saveIndexImp
( ostream& fileFile, 
  ostream& indexFile, 
  int& fileNumber )
{
  computeNumSequencesInFile(); // ensure have scanned to end of file
  fileFile << fileName_ << endl;
  SeqIndexInfo* pIndex = new SeqIndexInfo[seqPositions_.size()];
  for (int i(0) ; i < seqPositions_.size() ; i++)
  {
    pIndex[i].fileNum=fileNumber;
    pIndex[i].seqPos=seqPositions_[i];
  } // ~for

  indexFile.write( (const char*)pIndex, 
                   seqPositions_.size()*sizeof(SeqIndexInfo) );
  delete [] pIndex;
} // ~SequenceReaderFile::saveIndex



// End of file SequenceReaderFasta.cpp

---