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

// Description:
// Implements C interface for SSAHA


// Includes:
#include "SSAHAWrapper.h"
#include "SequenceReaderFasta.h"
#include "HashTableTranslated.h"
#include "SequenceEncoder.h"
#include <vector>
#include <string>

// ### Function Definitions ###

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

struct CustomHitStore : public vector<SSAHAHit>, public HitList
{
  virtual void addHit
   ( const PositionInDatabase& hitPos, 
      const SequenceOffset& queryPos )
  {
    push_back();
    back().seqNum = hitPos.sequence;
    back().seqPos = hitPos.offset;
  } // ~addHit

}; // ~class HitStore


class HashTablePackedCustom;
static HashTableGeneric* pHashTable(NULL);
static const int wordLength=5;
static SequenceEncoderProtein encoder(wordLength, cerr);
static CustomHitStore hits;
//vector<SSAHAHit> hits;
static PackedHitStore hitBuffer;
static vector<int> hitStarts; 
static bool isPacked=false;

class HashTablePackedCustom : public HashTablePackedProtein
{

 public:
  HashTablePackedCustom( ostream& monitoringStream=cerr, string name = ""):
    HashTablePackedProtein( monitoringStream, name ) 
    {}

  void convertHits( void );

};


// Name: HashTablePackedCustom::convertHits
// Based on HashTablePacked::convertHits
// Takes contents of hitBuffer, unpacks them, puts in hit store
void HashTablePackedCustom::convertHits( void )
 //( PackedHitStore& packedHits, HitList& hitListFwd )
{

  //  sort(hitBuffer.begin(),hitBuffer.end());

  vector<HitPacked>::iterator i(hitBuffer.begin());
  vector<SeqStartPos>::iterator j(seqStarts_.begin());
  //  hitListFwd.reserve(hitListFwd.size()+hitBuffer.size());
  hits.reserve( hits.size() + hitBuffer.size());

  // NB seqStarts_ must not be empty else call to back() segfaults

  while ( (j!=&seqStarts_.back()) && (i!=hitBuffer.end()) )
  {
    vector<SeqStartPos>::iterator ub(upper_bound(j,seqStarts_.end(),i->first));
    j=ub; j--;

    while ((i!=hitBuffer.end())&&(i->first<*ub))
    {
      //      cout << i->first << " ... " 
      //   << ub-seqStarts_.begin() << " " 
      //   << stepLength_*(i->first - *j) 
      //   << endl;
      //      hitListFwd.addHit( ub-seqStarts_.begin(), 
      //                    stepLength_*(i->first - *j),
      //                    i->second );
      hits.push_back();
      hits.back().seqNum= ub-seqStarts_.begin();
      hits.back().seqPos= stepLength_*(i->first - *j);
      //      hitListFwd.addHit( 
      //                    
      //                    i->second );
      
      i++;
      
    } // ~while
    j=ub;  
  } // ~while

} // ~HashTablePackedCustom::convertHits( void )




/* return name of the seqNum-th sequence in the table */
const char* getSubjectName( unsigned int seqNum )
{
  return pHashTable->getSequenceName( seqNum );
} // ~const char* getSubjectName( unsigned int seqNum )



int loadTable( const char* const tableName )
{
  clearSearch();
  clearTable();

  if (isPacked)
  {
    pHashTable = new HashTablePackedCustom( cerr, (string)tableName);
  }
  else
  {
    pHashTable = new HashTableFred( cerr, (string)tableName);
  }

  pHashTable->loadHashTable(); 
  // will throw exception if fails, eg if trying to load wrong table
  assert(pHashTable->getWordLength()==wordLength);
  pHashTable->setMaxNumHits( 100000000 );
}

/* perform a SSAHA search */
int doSearch
( 
 const char* const pWords,
 const int numWords,
 SSAHAHit** const pHits,
 int ** const pStarts
)
{
  clearSearch();
  //  SequenceEncoderProtein encoder(wordLength, cerr);
  WordSequence w;

  encoder.linkSeq(w);
  encoder.encode(pWords, numWords*wordLength);
  encoder.unlinkSeq();

  if(w.empty())
  {
    hitStarts.push_back(-1);
    *pHits=NULL;
    *pStarts=&hitStarts[0];
    return -1;
  }


  if (isPacked)
  {
    HashTablePackedCustom* pTable
      = dynamic_cast<HashTablePackedCustom*>(pHashTable);
    assert(pTable!=NULL);
    pTable->setQueryProtein();
    for (WordSequence::iterator i(w.begin());i!=&w.back();i++)
    {
      //    cout << "got: " << printResidue(*i, wordLength) << endl;
      hitBuffer.clear();
      pTable->matchWord( *i, hitBuffer );
      //    pHashTable->convertHits( hitBuffer, hits );
      pTable->convertHits();
      hitStarts.push_back(hits.size());
    }
  } // ~for
  else
  {
    HashTableFred* pTable
      = dynamic_cast<HashTableFred*>(pHashTable);
    assert(pTable!=NULL);
    //    int n(hits.size());
    for (WordSequence::iterator i(w.begin());i!=&w.back();i++)
    {
      pTable->matchWord( *i, hits );
    }
    hitStarts.push_back(hits.size());
  }


  //  hits.push_back();
  //  hits.back().seqNum=1103; hits.back().seqPos=-1103;
  //  hitStarts.push_back(1);
  hitStarts.push_back(-1);
  *pHits=&hits[0];
  *pStarts=&hitStarts[0];
  return 0;
}

/* Clear search results from RAM (automatically done by loadTable) */
int clearSearch( void )
{
  hits.clear();
  hitStarts.clear();
}

/* Clear SSAHA tables from RAM (automatically done by loadTable) */
int clearTable( void )
{
  delete pHashTable;
  pHashTable = NULL;
}



// End of file SSAHAWrapper.cpp

---