QueryManager/MatchStore.h

Go to the documentation of this file.

/*  Last edited: Apr  4 10:15 2002 (ac2) */

// #######################################################################

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

// Include guard:
#ifndef INCLUDED_MatchStore
#define INCLUDED_MatchStore

// Description:

// Includes:
#include "GlobalDefinitions.h"
#include <string>
#include <algorithm>
#include <deque>

// NB it is good practise for #include statements in header files to be
// replaced by forward declarations if at all possible

// ### Class Declarations ###


class HashTableGeneric;
class SourceReader;


// Match and subclasses

class Match
{
 public:
  virtual ~Match() {}

  virtual SequenceNumber getSubjectNum( void ) const=0;
  //  virtual string         getSubjectName( void ) const=0;
  virtual const char*         getSubjectName( void ) const=0;
  virtual SequenceOffset getSubjectStart( void ) const=0;
  virtual SequenceOffset getSubjectEnd( void ) const=0;

  virtual SequenceNumber getQueryNum( void ) const=0;
  virtual string         getQueryName( void ) const=0;
  virtual SequenceOffset getQueryStart( void ) const=0;
  virtual SequenceOffset getQueryEnd( void ) const=0;
  
  virtual int  getQuerySize( void ) const=0;
  virtual int  getNumBases(void ) const=0;
  virtual bool isQueryForward( void ) const=0;
  virtual bool isSubjectForward( void ) const=0;
  
  virtual void           print( void ) const=0;
 

}; // ~class Match


// Struct Name: MatchImp
// Description: This contains details of a gapped alignment between 
// the query sequence under analysis and a sequence from the subject database. 

class MatchStoreImp;

class MatchImp : public Match
{
 public:

  virtual ~MatchImp() {}

  virtual SequenceNumber getSubjectNum( void ) const 
    { return subjectNum_; }
  //  virtual inline string  getSubjectName( void ) const;
  virtual inline const char*  getSubjectName( void ) const;
  virtual SequenceOffset getSubjectStart( void ) const 
    { return subjectStart_; }
  virtual SequenceOffset getSubjectEnd( void ) const 
    { return subjectEnd_; }

  virtual inline SequenceNumber getQueryNum( void ) const;
  virtual inline string         getQueryName( void ) const;
  virtual SequenceOffset getQueryStart( void ) const 
    { return queryStart_; }
  virtual SequenceOffset getQueryEnd( void ) const 
    { return queryEnd_; }
  
  virtual inline int  getQuerySize( void ) const; 
  virtual int  getNumBases(void ) const 
    { return numBases_; }
  virtual bool isQueryForward( void ) const 
    { return isQueryForward_; }
  virtual bool isSubjectForward( void ) const 
    { return isSubjectForward_; }
  
  virtual void           print( void ) const {}
  //


  MatchImp(  MatchStoreImp* myStore,
             SequenceNumber subjectNum,
             SequenceOffset numBases,
             SequenceOffset queryStart,
             SequenceOffset queryEnd,
             SequenceOffset subjectStart,
             SequenceOffset subjectEnd,
             bool isQueryForward,
             bool isSubjectForward ):
    myStore_( myStore ),
    subjectNum_( subjectNum ),  
    numBases_( numBases ),  
    queryStart_( queryStart ),
    queryEnd_( queryEnd ),
    subjectStart_( subjectStart ),
    subjectEnd_( subjectEnd ),
    isQueryForward_( isQueryForward ),
    isSubjectForward_( isSubjectForward ){}

  SequenceNumber subjectNum_;
  SequenceOffset numBases_;

  SequenceOffset queryStart_;
  SequenceOffset queryEnd_;
  SequenceOffset subjectStart_;
  SequenceOffset subjectEnd_;

  bool isQueryForward_;
  bool isSubjectForward_;

  MatchStoreImp* myStore_;

};


// MatchStore and subclasses

class MatchStore : public vector<Match*>
{
 public:
  virtual ~MatchStore() {}

  virtual void clear( void ) {}


  virtual void printResult( ostream& outputStream = cout ) const {}

  virtual void addMatch(  SequenceNumber subjectNum,
                          SequenceOffset numBases,
                          SequenceOffset queryStart,
                          SequenceOffset queryEnd,
                          SequenceOffset subjectStart,
                          SequenceOffset subjectEnd,
                          bool isQueryForward,
                          bool isSubjectForward ) {}
};

class HitListVector;
class NameReader;

class MatchStoreImp : public MatchStore
{
  friend class MatchImp;
  public:
  MatchStoreImp( string queryName,
                 int        queryNum,
                 int        queryBases,
                 const NameReader& reader) :
  queryName_ ( queryName ),
  queryNum_( queryNum ),
  queryBases_( queryBases ), 
  reader_( reader )
    {
      //      reserve(1000);
      //     imps_.reserve(1000);
    }

  ~MatchStoreImp();
  
  virtual void addMatch(  SequenceNumber subjectNum,
                          SequenceOffset numBases,
                          SequenceOffset queryStart,
                          SequenceOffset queryEnd,
                          SequenceOffset subjectStart,
                          SequenceOffset subjectEnd,
                          bool isQueryForward,
                          bool isSubjectForward );

  // Function Name: printResult
  // Arguments: const QueryResult& (in), ostream& (out)
  // Convert the contents of result to ASCII and send to the specified
  // output stream
  void printResult
  ( ostream& outputStream = cout ) const;

  virtual void clear( void ); 

  void setup( void );

 protected:
  vector<MatchImp*> imps_;

  string     queryName_;
  int        queryNum_;
  int        queryBases_;
  //  HashTable& subjectTable_;
  //  mutable map< SequenceNumber, string > names_;
  const NameReader& reader_;

}; // ~MatchStore



// MatchTask and subclasses

class MatchTask
{
 public:
  virtual ~MatchTask() {}
  virtual void operator()( MatchStore& store ) =0;
};

template <class TASK1, class TASK2> class CombineTask : public MatchTask
{
  public:
    CombineTask(TASK1& t1, TASK2& t2) : t1_(t1), t2_(t2) {}
    CombineTask( void ) : t1_imp(), t2_imp(), t1_(t1_imp), t2_(t2_imp) {}
    virtual void operator()( MatchStore& store )
{
    t1_(store);
    t2_(store);
  }
  private:
  TASK1 t1_imp;
  TASK2 t2_imp;
  TASK1& t1_;
  TASK2& t2_;
};

// CombineTaskVirtual uses virtual function polymorphism instead of templates:
// more versatile as you don't have to know exact type of task1 and
// task2 at compile time. Cost of virtual function likely to be negligible
// wrt to tasks themselves (sorting, printing etc)

class CombineTaskVirtual : public MatchTask
{
 public:
  CombineTaskVirtual(MatchTask& t1, MatchTask& t2) : t1_(t1), t2_(t2) {}
  virtual void operator()( MatchStore& store )
  {
    t1_(store);
    t2_(store);
  }
  private:
  MatchTask& t1_;
  MatchTask& t2_;
};

// prints out all matches
class MatchTaskPrint : public MatchTask
{
 public: 
  MatchTaskPrint( ostream& outputStream = cout ) :
    outputStream_( outputStream ) {}
  virtual void operator()( MatchStore& store );
 private:
  ostream& outputStream_;

};

// prints out all matches in `parser friendly format'
class MatchTaskPrintTabbed : public MatchTask
{
 public: 
  MatchTaskPrintTabbed( ostream& outputStream = cout ) :
    outputStream_( outputStream ) {}
  virtual void operator()( MatchStore& store );
 private:
  ostream& outputStream_;
};

// same as MatchTaskPrint, except reverses frame of query
// coordinates for reverse matches
class MatchTaskPrintReverse : public MatchTask
{
 public: 
  MatchTaskPrintReverse( ostream& outputStream = cout ) :
    outputStream_( outputStream ) {}
  virtual void operator()( MatchStore& store );
 private:
  ostream& outputStream_;

};

// same as MatchTaskPrintTabbed, except reverses frame of query
// coordinates for reverse matches
class MatchTaskPrintTabbedReverse : public MatchTask
{
 public: 
  MatchTaskPrintTabbedReverse( ostream& outputStream = cout ) :
    outputStream_( outputStream ) {}
  virtual void operator()( MatchStore& store );
 private:
  ostream& outputStream_;
};


class SortByMatchSize 
{
public:
  SortByMatchSize() {}
  bool operator()( const Match* p1, const Match* p2 ) const
  { 
    // We want the largest match at the top
    return ( (p1->getNumBases()) > (p2->getNumBases() ) );
  }
};

class SortBySubjectName 
{
 public:
  SortBySubjectName() {}
  bool operator()( const Match* p1, const Match* p2 ) const
  { 
    return ( (p1->getSubjectName()) < (p2->getSubjectName() ) );
  }
};

class SortByQueryStart
{
 public:
  SortByQueryStart() {}
  bool operator()( const Match* p1, const Match* p2 ) const
  { 
    return ( (p1->getQueryStart()) < (p2->getQueryStart() ) );
  }
};

class SortByDiff
{
 public:
  SortByDiff() {}
  bool operator()( const Match* p1, const Match* p2 ) const
  { 
    return (   (p1->getQueryStart()-p1->getSubjectStart())
             > (p2->getQueryStart()-p2->getSubjectStart())  );
  }
};

class SortByPercentMatch
{
 public:
  SortByPercentMatch() {}
  bool operator()( const Match* p1, const Match* p2 ) const
  { 
    return (   (p1->getNumBases()*(p2->getQueryEnd()-p2->getQueryStart()+1))
             > (p2->getNumBases()*(p1->getQueryEnd()-p1->getQueryStart()+1)));
  }
};
  
class SortNull
{
public:
  SortNull() {}
  bool operator()( const Match* p1, const Match* p2) const 
  {
    return false;
  }
};

#ifdef OLD_SORTBYMATCH_DEFINITION
class SortByMatch 
{
public:
  SortByMatch() {}
  bool operator()( const Match* p1, const Match* p2 ) const
  { 
    return (    (p1->getNumBases()) > (p2->getNumBases() )
             || (    (p1->getNumBases()) == (p2->getNumBases() )
                  && (p1->getSubjectName()) < (p2->getSubjectName() ) )
             || (    (p1->getSubjectName()) == (p2->getSubjectName() )
                  && (p1->getQueryStart()) == (p2->getQueryStart() ) ) );
  }
};
#endif

// Sorting by subject name is very slow, so sort by subject number instead
// TC 25.3.2
class SortByMatch 
{
public:
  SortByMatch() {}
  bool operator()( const Match* p1, const Match* p2 ) const
  { 
    return (    (p1->getNumBases()) > (p2->getNumBases() )
             || (    (p1->getNumBases()) == (p2->getNumBases() )
                  && (p1->getSubjectNum()) < (p2->getSubjectNum() ) )
             || (    (p1->getSubjectNum()) == (p2->getSubjectNum() )
                  && (p1->getQueryStart()) == (p2->getQueryStart() ) ) );
  }
};


template < class C1, class C2, class C3=SortNull > 
class CombineSort
{
public:
  CombineSort() : c1_(), c2_(), c3_() {}
  bool operator()( const Match* p1, const Match* p2) const
  {
    return (    c1_(p1,p2)
             || ( (!c1_(p2,p1)) && c2_(p1,p2) ) 
             || ( (!c2_(p2,p1)) && c3_(p1,p2) ) );
  }
  C1 c1_;
  C2 c2_;
  C3 c3_;
};

template <class SORTER> class MatchTaskSort : public MatchTask
{
public:
  MatchTaskSort
  ( 
    unsigned int maxToSort = 1<<30,
    double partialThreshold = 0.0
  ) : sorter_(), maxToSort_(maxToSort), partialThreshold_(partialThreshold) {}
  // TBD should be virtual???
  void operator()(MatchStore& store )
  {
    if (    ( store.size() <= maxToSort_)
         || (  ((double)maxToSort_) / ((double)store.size())
               > partialThreshold_  )  )
      sort( store.begin(), store.end(), sorter_ );
    else 
    {
      partial_sort( store.begin(), 
                    store.begin()+maxToSort_, 
                    store.end(), 
                    sorter_ );
    }
    store.resize(min((unsigned int)store.size(),maxToSort_));
  }
private:
  unsigned int maxToSort_;
  double partialThreshold_;
  SORTER sorter_;
}; // class MatchSort

#ifdef MOVED_TO_ALIGNERDOTH
// Declarations for MatchTaskAlign

typedef int ScoreType;
typedef ScoreType ScoreTable[65536];
typedef unsigned char uchar;
typedef unsigned char PathType;

typedef unsigned char AlignBreakType;

// alignment of two sequences is described by a vector of AlignInfo 
// Meant for sequences that contain large stretches of matches
// punctuated by the occasional AlignBreak

struct AlignInfo
{
  AlignInfo( void ) : numMatches(0), breakType(0) {}
  int numMatches;
  AlignBreakType breakType;
};


//void fillScoreTable( ScoreTable& table );

class MatchTaskAlign : public MatchTask
{
public:
  MatchTaskAlign( SourceReader& querySource, 
                  SourceReader& subjectSource,
                  int numCols,
                  ostream& outputStream,
                  bool reverseQueryCoords, 
                  bool doAlignment );
  virtual ~MatchTaskAlign();

  virtual void operator()(MatchStore& store);

  // TBD rest of these could be made private

  void align( const char* pQuery, int queryStart, int queryEnd,
              const char* pSubject, int subjectStart, int subjectEnd );
  void createAlignment
  ( deque<AlignInfo>& alignment,
    const char* pQuery, int queryStart, int queryEnd,
    const char* pSubject, int subjectStart, int subjectEnd );
  void formatAlignment
  ( deque<AlignInfo>& alignment,
    const char* pQuery, int queryStart, int queryEnd,
    const char* pSubject, int subjectStart, int subjectEnd );

  bool matchChar( const char a, const char b )
  {
    return (tolower(a)==tolower(b));
  }

  virtual void outputAlignmentLine( void );

  void reverseQueryCoords( bool yesOrNo= true) 
  { 
    reverseQueryCoords_ = yesOrNo;
  }

  void doAlignment( bool yesOrNo= true) 
  { 
    doAlignment_ = yesOrNo;
  }
 
protected:

  ostream& outputStream_;
  SourceReader& querySource_;
  SourceReader& subjectSource_;
  // reverseQueryCoords_ - if `true' reverses orientation of
  // query coordinates, so they are as seen in the forward frame
  // Basically set to true for ssaha executable but false in 
  // ssahaClient, because in that case any necessary reversals
  // are done in ssahaServer
  bool reverseQueryCoords_;
  // doAlignment_ - if `true' prints match description line plus
  // graphical alignment for each match, if `false' just prints
  // match description line
  bool doAlignment_;

  char* pBufSeq1;
  char* pBufAlign;
  char* pBufSeq2;
  int numCols_;

  ScoreTable* pTable_;


};

// Class Name: MatchTaskAlignDNA
// Description: Aligns using DNA score table
class MatchTaskAlignDNA : public MatchTaskAlign
{
 public:
  MatchTaskAlignDNA( SourceReader& querySource, 
                     SourceReader& subjectSource,
                     int numCols=80,
                     ostream& outputStream=cout,
                     bool reverseQueryCoords=true,
                     bool doAlignment=true);
 protected:
  void fillScoreTable( void );

  static bool isTableFilled_;
  static ScoreTable table_;
};

//#ifdef DATA_SOURCES_FOR_BLOSUM62_MATRIX
The below is cut and pasted from
http://www.ncbi.nlm.nih.gov/IEB/ToolBox/C_DOC/lxr/source/data/BLOSUM62

-- quote --

  1 #  Matrix made by matblas from blosum62.iij
  2 #  * column uses minimum score
  3 #  BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
  4 #  Blocks Database = /data/blocks_5.0/blocks.dat
  5 #  Cluster Percentage: >= 62
  6 #  Entropy =   0.6979, Expected =  -0.5209
  7    A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  *
  8 A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1  0 -4 
  9 R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1  0 -1 -4 
 10 N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  3  0 -1 -4 
 11 D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4  1 -1 -4 
 12 C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4 
 13 Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0  3 -1 -4 
 14 E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4 
 15 G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -2 -1 -4 
 16 H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0  0 -1 -4 
 17 I -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3 -3 -3 -1 -4 
 18 L -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1 -4 -3 -1 -4 
 19 K -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2  0  1 -1 -4 
 20 M -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1 -3 -1 -1 -4 
 21 F -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1 -3 -3 -1 -4 
 22 P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2 -2 -1 -2 -4 
 23 S  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2  0  0  0 -4 
 24 T  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0 -1 -1  0 -4 
 25 W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3 -4 -3 -2 -4 
 26 Y -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1 -3 -2 -1 -4 
 27 V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4 -3 -2 -1 -4 
 28 B -2 -1  3  4 -3  0  1 -1  0 -3 -4  0 -3 -3 -2  0 -1 -4 -3 -3  4  1 -1 -4 
 29 Z -1  0  0  1 -3  3  4 -2  0 -3 -3  1 -1 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4 
 30 X  0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2  0  0 -2 -1 -1 -1 -1 -1 -4 
 31 * -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4  1 

-- unquote --

In addition:
Score of any character with - is -10
Score of - with - is -1
These are the gap open and gap extension penalties used with BLOSUM62 in
Huang/Zhang, CABIOS 12(6) 1996, pp. 497-506

//#endif

static const char blosumResidues[] 
= "ARNDCQEGHILKMFPSTWYVBZX";
// 12345678901234567890123

static const ScoreType blosumScores[] =
{
4, // A
-1, 5, // R
-2, 0, 6, // N
-2, -2, 1, 6, // D
0, -3, -3, -3, 9,// C
-1, 1, 0, 0, -3, 5,// Q
-1, 0, 0, 2, -4, 2, 5,// E
0, -2, 0, -1, -3, -2, -2, 6, // G
-2, 0, 1, -1, -3, 0, 0, -2, 8, // H
-1, -3, -3, -3, -1, -3, -3, -4, -3, 4, // I
-1, -2, -3, -4, -1, -2, -3, -4, -3, 2, 4, // L
-1, 2, 0, -1, -3, 1, 1, -2, -1, -3, -2, 5, // K
-1, -1, -2, -3, -1, 0, -2, -3, -2, 1, 2, -1, 5, // M
-2, -3, -3, -3, -2, -3, -3, -3, -1, 0, 0, -3, 0, 6, // F
-1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4, 7, // P
1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2, 0, -1, -2, -1, 4, // S
0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 1, 5, // T
-3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1, 1, -4, -3, -2, 11, // W
-2, -2, -2, -3, -2, -1, -2, -3, 2, -1, -1, -2, -1, 3, -3, -2, -2, 2, 7,// Y
0, -3, -3, -3, -1, -2, -2, -3, -3, 3, 1, -2, 1, -1, -2, -2, 0, -3, -1, 4,// V
-2, -1, 3, 4, -3, 0, 1, -1, 0, -3, -4, 0, -3, -3, -2, 0, -1, -4, -3, -3, 4,// B
-1, 0, 0, 1, -3, 3, 4, -2, 0, -3, -3, 1, -1, -3, -1, 0, -1, -3, -2, -2, 1, 4, // Z
0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, 0, 0, -2, -1, -1, -1, -1, -1 //X
};

static const int numBlosumResidues = 23;

// Class Name: MatchTaskAlignProtein
// Description: Aligns protein using score table based on BLOSUM62
class MatchTaskAlignProtein : public MatchTaskAlign
{
 public:
  MatchTaskAlignProtein( SourceReader& querySource, 
                         SourceReader& subjectSource,
                         int numCols=80,
                         ostream& outputStream=cout );
 protected:
  void fillScoreTable( void );

  static bool isTableFilled_;
  static ScoreTable table_;
};

#endif







// Class Name :
// Description: 

// ### Function Declarations ###

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

// End of include guard:
#endif

// End of file MatchStore.h

---