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

// Description:

// Includes:

#include "SequenceReaderFilter.h"
#include <strstream>

// ### Function Definitions ###

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

  // Function Name: Constructor
  // Arguments: ostream&
  // Takes the data from seqFile and places it in seqData_, seqBasesInLast_
  // and seqNames_

SequenceReaderFilter::SequenceReaderFilter
( SequenceReader* pSeq, 
  ifstream* pFilterSource,
  ostream& monitoringStream ) :
  SequenceReader( monitoringStream ),
  pSeq_(pSeq),
  pFilterNames_( new StringHash ),
  numFiltered_(0)
{
  monitoringStream_ 
    << "constructing SequenceReaderFilter" << endl;
  readFilterNames( pFilterSource );
} // ~constructor

SequenceReaderFilter::SequenceReaderFilter
( SequenceReader* pSeq, 
  const char* filterFileName, 
  ostream& monitoringStream ) :
  SequenceReader( monitoringStream ),
  pSeq_(pSeq),
  pFilterNames_( new StringHash ),
  numFiltered_(0)
{
  monitoringStream_ 
    << "constructing SequenceReaderFilter for file " 
    << filterFileName << endl;
  readFilterNames( new ifstream( filterFileName ) );
} // ~constructor

// Function Name: readFilterNames
// Arguments: ifstream*
// Returns: void
// Reads in a list of names to be filtered from an ifstream*
// into set<string>*. Once done, deletes the ifstream
void SequenceReaderFilter::readFilterNames( ifstream* pFilterSource )
{
  if ( pFilterSource->fail() )
  {
    throw SSAHAException
      ( "SequenceReaderFile - unable to open filter file ");

  } // ~if

  //  string thisName;

  //  while (*pFilterSource>>thisName) pFilterNames_->insert(thisName);
  pFilterNames_->push_back( (string) "" );
  while (*pFilterSource>>pFilterNames_->back()) pFilterNames_->push_back( (string) "" );
  pFilterNames_->pop_back();
  pFilterNames_->makeBins();

  monitoringStream_ << "Read in " << pFilterNames_->size() 
                    << " names to filter out" << endl;
  filterNums_.push_back(0);
  pFilterSource->close();
  delete pFilterSource;
}




// Function Name: Copy constructor
// Arguments:
SequenceReaderFilter::SequenceReaderFilter( const SequenceReaderFilter& rhs ):
SequenceReader( rhs.monitoringStream_ ),
pSeq_( rhs.pSeq_->clone() )
{
  monitoringStream_ << "copy constructing SequenceReaderFilter" << endl;
  // copy set across
  monitoringStream_ << "copy constructor not implemented!!" << endl;
  assert(1==0);
} // ~destructor

  // Function Name: Destructor
  // Arguments:
SequenceReaderFilter::~SequenceReaderFilter()
{
  monitoringStream_ << "destructing SequenceReaderFilter" << endl;
  delete pSeq_;
  delete pFilterNames_;
} // ~destructor

bool SequenceReaderFilter::findSequence( SequenceNumber seqNum )
{
  if ( seqNum <= filterNums_.size()-1 ) return true; // no neeed!
  if (allSequencesRead_)
  {
    monitoringStream_ 
      << "Requested sequence number (" << seqNum 
      << ")\nis outside range of sequences in file (1 to " 
      << numSequencesInFile_ << ")." << endl;
    return false;
  } // ~if

  SequenceNumber lastInSource(filterNums_.back());
  // = (filterNums_.empty() ? 0 : filterNums_.back()); no need, never empty
  std::ostrstream buffer;
  string name;

  while (filterNums_.size()!=seqNum+1)
  {

    ++lastInSource;

    if ( pSeq_->printName( buffer, lastInSource ) == false )
    {
      if(allSequencesRead_==false)
      {
        allSequencesRead_=true;
        if (pFilterNames_!=NULL) 
        {
          monitoringStream_ 
          << "deleting filter name data, no longer needed" << endl;
          delete pFilterNames_;
          pFilterNames_=NULL;  
        } // ~if
        numSequencesInFile_ = filterNums_.size()-1;
        monitoringStream_ 
          << "Requested sequence number (" << seqNum
          << ")\nis outside range of sequences in file (1 to " 
          << numSequencesInFile_ << ")." << endl;

      } // ~if
      return false;
    } // ~catch

    name=buffer.str();
    string::size_type nameEnd=name.find_first_of('\t');
    // NB stop at tab not spc
    if (nameEnd!=string::npos) name.erase(nameEnd,name.size());

    //    if ( pFilterNames_->find(name)==pFilterNames_->end() )
    if ( !pFilterNames_->isPresent(name) )
    {
      filterNums_.push_back(lastInSource);
      //  monitoringStream_ << filterNums_.size()-1 << " " << lastInSource << endl;
    } // ~if
    else
    {
      //      monitoringStream_ << filterNums_.size()-1 << " " << lastInSource 
      //                << " Filtered: " << name << " " 
      //                << *(pFilterNames_->find(name)) << endl;
      numFiltered_++;
    } // ~else


    buffer.freeze(false); // release memory to write further stuff to buffer  
    buffer.seekp(0,ios::beg); // ensure next write is to start of buffer
  
  } // ~while

  return true;

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

// 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 SequenceReaderFilter::getNextSequence( WordSequence& nextSeq, 
                                           int wordLength )
{
  if (!findSequence(lastSequenceNumber_+1)) 
  {
    assert(filterNums_.size()==lastSequenceNumber_+1);
    return -1;
  } // ~if
  else 
  {
    ++lastSequenceNumber_;
    return pSeq_->getSequence( nextSeq, filterNums_[lastSequenceNumber_], 
                                 wordLength );
  } // ~else
} // ~SequenceReaderFilter::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 SequenceReaderFilter::getSequence
( WordSequence& nextSeq, SequenceNumber sequenceNumber, int wordLength )
{
  if (!findSequence(sequenceNumber)) return -1;
  assert(filterNums_.size()>=sequenceNumber+1);
  lastSequenceNumber_=sequenceNumber;
  return pSeq_->getSequence
    ( nextSeq, filterNums_[sequenceNumber], wordLength );
  
} // ~SequenceReaderFilter::getSequence


void SequenceReaderFilter::rewind( void ) 
{
  pSeq_->rewind();
  lastSequenceNumber_ = 0;
} // ~rewind 

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

// Function Name: getSequenceName
// Arguments: string& (out), SequenceNumber (in)
// Returns:   void
// Fills a string with the name of the requested sequence
bool SequenceReaderFilter::printName( ostream& os, SequenceNumber seqNum) 
{
  if (!findSequence(seqNum)) return false;
  assert(filterNums_.size()>=seqNum+1);
  pSeq_->printName(os, filterNums_[seqNum]);
  return true;
} // ~SequenceReaderFilter::printName( ostream& os, SequenceNumber seqNum)  

// Function Name: printSideInfo
// Arguments: string& (out), SequenceNumber (in)
// Returns:   void
// Fills a string with the name of the requested sequence
bool SequenceReaderFilter::printSideInfo
( ostream& os, SequenceNumber seqNum )
{
  if (!findSequence(seqNum)) return false;
  assert(filterNums_.size()>=seqNum+1);
  pSeq_->printName(os, filterNums_[seqNum]);
  return true;
} // ~bool SequenceReaderFilter::printSideInfo

// Function Name: printSource
// Arguments: string& (out), SequenceNumber (in)
// Returns:   void
// Fills a string with the name of the requested sequence
bool SequenceReaderFilter::printSource
( ostream& os, SequenceNumber seqNum )
{
  if (!findSequence(seqNum)) return false;
  assert(filterNums_.size()>=seqNum+1);
  pSeq_->printName(os, filterNums_[seqNum]);
  return true;
} // ~SequenceReaderFilter::printSource

SequenceNumber SequenceReaderFilter::computeNumSequencesInFile( void )
{
  try
  {
    // this will always work out num seqs in pSeq_ 
    findSequence( 1+pSeq_->getNumSequencesInFile() );
  } // ~try
  catch( const NumberOutOfRange& err )
    {} // ~catch

  return numSequencesInFile_;

} // ~SequenceReaderFilter::computeNumSequencesInFile( void )

// Functions to implement SourceReader functionality

// Function Name: extractSource
// This extracts the source data for bases seqStart to seqEnd inclusive
// of sequence seqNum and places it in source
void SequenceReaderFilter::extractSource
( char** pSource, 
  SequenceNumber seqNum,
  SequenceOffset seqStart,
  SequenceOffset seqEnd )
{
  //  cout << seqNum << " " << filterNums_[seqNum] << " " << seqStart << " " << seqEnd << endl; 
  pSeq_->extractSource( pSource, filterNums_[seqNum], seqStart, seqEnd );
} // ~void SequenceReaderFilter::extractSource


// Function Name: saveIndexImp
// This makes and saves the index for *this, by kinda monkeying around
// with the index for *pSeq_ 
void SequenceReaderFilter::saveIndexImp
( ostream& fileFile, 
  ostream& indexFile, 
  int& fileNumber )
{
  std::strstream indexTemp;
  pSeq_->saveIndexImp( fileFile, indexTemp, fileNumber );
  const SeqIndexInfo* pIndex= (const SeqIndexInfo*)indexTemp.str();

  // writing individually is slow as a slow thing on Valium ...
  //  int numSeqs(getNumSequencesInFile());
  //  for (int i(1); i<=numSeqs; i++)
  //  {
  //    //        cout << "Index: " << i << " " << filterNums_[i]-1 << " "  
  //    //     << pIndex[filterNums_[i]-1].fileNum
    //   << " " << pIndex[filterNums_[i]-1].seqPos << endl;
  //   indexFile.write
  //    ( (const char*)&pIndex[ filterNums_[i]-1 ], sizeof(SeqIndexInfo) );
  //  } // ~for

  int numSeqs(getNumSequencesInFile());
  vector<SeqIndexInfo> newIndex;
  newIndex.reserve(numSeqs);

  for (int i(1); i<=numSeqs; i++)
  {
    //        cout << "Index: " << i << " " << filterNums_[i]-1 << " "  
    //     << pIndex[filterNums_[i]-1].fileNum
    //   << " " << pIndex[filterNums_[i]-1].seqPos << endl;
    newIndex.push_back(pIndex[ filterNums_[i]-1 ]);
  } // ~for

  indexFile.write
    ( (const char*)&newIndex[ 0 ], numSeqs*sizeof(SeqIndexInfo) );



} // ~void SequenceReaderFilter::saveIndexImp


// StringHash member functions


void StringHash::makeBins( int numBins )
{
  //  cout << "Hashing " << size() << " sequence names into "
  //   << numBins << " bins" << endl;
  
  numBins_ = numBins;
    bins_.clear();
    bins_.resize(numBins);
    for (vector<string>::const_iterator i(begin());i!=end();i++)
    {
      (bins_[ H_(i->c_str()) % numBins ]).push_back((string*)&(*i));
      //      (bins_[ H_(i->c_str()) % numBins ]).push_back();
      //  (bins_[ H_(i->c_str()) % numBins ]).back() = i;
    }
} // StringHash

size_t StringHash::makeNumBins( void ) const
{
  static const size_t primes[] = { 103,1009,7013,10007,88883,100043 };
  static const int numPrimes = 6;
    size_t minSize = size() + (size()>>3); // current size +12.5%
    for (int i(0);i<numPrimes;i++)
      if (minSize<primes[i]) return primes[i];
    minSize |=1; // make sure it's odd
    if (minSize%5==0) minSize+=2; // make sure it's not a mult of 5
    // (cos then hash fn gives lousy performance)
    return minSize;
}






// End of file SequenceReaderFilter.cpp

---