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

// Description:

// Includes:

#include "GlobalDefinitions.h"
#include "HashTable.h"
#include "HashTablePacked.h"
#include "HashTableTranslated.h"
#include "SequenceReaderString.h"
#include "SequenceReaderFasta.h"
#include "SequenceEncoder.h"
#include "GenerateTestFastaFiles.h"
#include "MatchStoreUngapped.h"
#include "MatchStore.h"
#include "QueryManager.h"
#include <assert.h>
#include <string>

// ### Function Definitions ###

// breaks encapsulation for debugging/test purposes
class HashTableTranslatedPeek : public HashTableTranslated
{
public:
  HashTableTranslatedPeek( void ): HashTableTranslated(cout) {}
  HashTableComponent& peekFwd( void ) { return hashFwd_; }
  HashTableComponent& peekRev( void ) { return hashRev_; }
  
};




int main( void )
{

  try {

  int numTests = 0;

  cout << "*******************************************" << endl << endl;
  cout << "     Test of class HashTable" << endl << endl;
  cout << "*******************************************" << endl << endl;

  int numSeqs = 10;
  int seqSize = 100;
  int wordLength = 10;
  int maxHits = 50;

  // Generate a random sequence of (numSeqs*seqSize) base pairs ...
  // 1128 is the seed value for the random number generator
  BaseGenerator testBases( numSeqs * seqSize, 1128 );

  // ... make a fasta file from this sequence, breaking the sequence into
  // sequences of seqSize base pairs, with no overlaps (i.e. numSeqs of them)
  // Haven't specified a file name, so will default to test_subject.fasta 
  testBases.generateSubjectFile(seqSize,0);

  SequenceReaderFasta testReader("test_subject.fasta",cout);

  // ---

  cout << "Test " << ++numTests <<": test of function getNextSequence" 
       << endl << endl;

  HashTable  testHash(cout,"test_save");
  HashTableFactory creator(cout);


  // check init flag is not set on construction ...
  assert(testHash.isInitialized() == false );

  creator.createHashTable(testHash,testReader,wordLength,maxHits);

  // ... but is set once the hash table has been set up
  assert(testHash.isInitialized() == true );

  testReader.rewind();

  HitListVector hits;

  WordSequence  seq;

  string name1, name2, name3;

  // The next bit takes sequence data from testReader and uses it to
  // query the hash table. The hash table was also created from testReader,
  // so we are checking that all sequence data 'finds itself' in the hash
  // table in the correct position

  for ( int i(1) ; i <= numSeqs ; i++ )
  { // for each sequence in testReader ...
    testReader.getNextSequence(seq,wordLength);

    // ... check that the name strings match 
    testReader.getLastSequenceName(name1);
    testHash.getSequenceName(name2,i);
    //    assert(testHash.getSequenceName(i)==
    //   testHash.getNameReader().getSequenceName(i));
    name3 = (string) testHash.getSequenceName(i);
    assert(name1 == name2 );
    assert(name2 == name3 );
    cout << seq.getNumBasesInLast() << "!!\n";
    // ... go through the Words in the sequence one by one and look
    //    for matches in the hash table
    for ( int j(0) ; j < (seqSize/wordLength) ; j++ )
    {
       cout << j << " " << printBase(seq[j],wordLength) << endl; 
      
       assert((seq[j]&gCursedWord)==0);
       testHash.matchWord(seq[j],hits,0);      
       // ... check we have exactly one hit for the word ...
       cout << hits.size() <<endl;
       assert(hits.size() == 1);
       //       cout << i << " " << j << " " 
       //          << hits[0].sequence << " " 
       //         << hits[0].offset << endl; 
       // ... and that the sequence number and offset for the hit are OK

       assert(hits[0].subjectNum == i );
       assert(hits[0].diff == j*wordLength );
       assert(hits[0].queryPos == 1);
 
       hits.clear();
 
    } // ~for j

    // Now try matchWord for the sequence as a whole
    testHash.matchWord(seq,hits); 
    cout << hits.size() << "!!!!!\n";
    assert( hits.size() == (seqSize/wordLength) );
    for ( int j(0) ; j < hits.size() ; j ++ )
    {

      assert(hits[j].subjectNum == i);
      assert(hits[j].diff   == 0);
      assert(hits[j].queryPos == 1 + ( j * wordLength ) );
     
    } // ~for j

    hits.clear();
    seq.clear();
  } // ~for i


  cout << "Test passed!" << endl << endl; 

  // ---

  cout << "Test " << ++numTests <<": test of save and load of hash table" 
       << endl << endl;

  // save hash table to files test_save.body, test_save.head, test_save.name
  //  testHash.saveHashTable("test_save");
  //  creator.saveHashTable(testHash,"test_save");
  creator.saveHashTable(testHash);

  // Create another hash table ...
  HashTable testLoad(cout,"test_save");

  // ... and initialize it from these saved files
  //  testLoad.loadHashTable("test_save");
  creator.loadHashTable(testLoad);

  // Check init flag has been set ...
  assert( testLoad.isInitialized() == true );
  
  // and that word length matches original
  assert( testLoad.getWordLength() == wordLength );

  testLoad.setMaxNumHits( testHash.getMaxNumHits() );

  assert( testHash.getMaxNumHits()==testLoad.getMaxNumHits() );

  HitListVector hitsOrig, hitsCopy;

  testReader.rewind();

  // Now we run all the sequences in testReader as queries on the old
  // and new hash tables, storing the resulting hits in hitsOrig and
  // hitsCopy respectively.

  for ( int i(0) ; i < numSeqs ; i++ )
  {

    // ... check that the name strings match 
    testHash.getSequenceName(name1,i+1);
    testLoad.getSequenceName(name2,i+1);
    assert(name1 == name2 );
    assert(name1 == (string)testHash.getSequenceName(i+1));
    assert(name1 == (string)testLoad.getSequenceName(i+1));
    // ... check that the sequence sizes are the same

    //    cout << i << " " << testHash.getSequenceSize(i+1) << " " 
    // << testLoad.getSequenceSize(i+1) << endl; // %%%%
    assert( testHash.getSequenceSize(i+1) == testLoad.getSequenceSize(i+1) );

    

    // ... add hit information to hitsOrig and hitsCopy
    seq.clear();
    testReader.getNextSequence(seq,wordLength);
    testHash.matchWord(seq,hitsOrig);
    testLoad.matchWord(seq,hitsCopy);
  } // ~for i

  


  // Note that neither testCopy or testOrig are cleared between queries,
  // so by this point they contain all the hit positions obtained for
  // all queries. We check that this info is identical for both sets.
  assert(hitsOrig.size() == hitsCopy.size() );
  assert(hitsOrig == hitsCopy );

  cout << "Test passed!" << endl << endl; 

  // ---

  cout << "Test " << ++numTests <<": further test of function matchWord" 
       << endl << endl;

  // Next few lines hash string testSeq into a hash table
  string testSeq = 
  "AGCTGGCTAGTCTAACGTGTCTGTAGATCGTACGTGCATGCATGGGCAGTCTCTGTATCACTGGGTGGCC";
// 1234567890123456789012345678901234567890123456789012345678901234567890
   
  SequenceReaderString shiftReader(testSeq,cout);

  HashTable   shiftHash(cout);
  creator.createHashTable(shiftHash,shiftReader,wordLength,maxHits);

  for ( int i(0) ; i < wordLength ; i++ )
  {
    testSeq = testSeq.substr(1); // delete first character
    { // braces ensure a new instance is created each time round loop
      SequenceReaderString queryReader(testSeq,cout);
      seq.clear();
      hits.clear();
      queryReader.getNextSequence(seq,wordLength);
      shiftHash.matchWord(seq,hits);

      if (i==(wordLength-1))
        assert( hits.size() == seq.size() - 1 );
        else assert( hits.size() == 0 );
    }
  } // ~for i

  cout << "Test passed!" << endl << endl; 

  // ---

  cout << "Test " << ++numTests 
       <<": test that only limited number of hits are stored" 
       << endl << endl;

  string tenAs = "AAAAAAAAAA";
  testSeq = "";
  for ( int i(0) ; i < 50 ; i++ ) testSeq += tenAs;

  // testSeq now consists of 500 'A' characters. Next we hash this
  // sequence. With a word length of 10 it consists of 50 words
  // 'AAAAAAAAAA', but because we set maxNumHits to 20, only 20 should
  // be stored 
  SequenceReaderString reader(testSeq,cout);
  HashTable hash(cout);
  creator.createHashTable(hash,reader,10,20);
  assert(hash.getMaxNumHits()==20);

  hits.clear();
  // Word 'AAAAAAAAAA' is encoded as zero in our coding scheme
  Word allAs(0);
  // matchWord should only place 20 lots of hit data into hits ...
  hash.matchWord(allAs,hits);
  // ... check this
  assert(hits.size()==0);
  cout << hash.size(allAs) << endl;
  //  assert(hash.size(allAs)==50);
  hash.setMaxNumHits(50);
  assert(hash.getMaxNumHits()==50);
  hash.matchWord(allAs,hits);
  assert(hits.size()==50);

  cout << "Test passed!" << endl << endl; 

  // ---

  cout << "Test " << ++numTests 
       <<": test of tandem repeat masking\n\n"; 

  // Each sequence string in repeat has 90 bases. ith sequence repeats
  // every i bases.
  string repeats[] = 
  {
  "GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG",
  "GCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCGCG",
  "AGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA",
  "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACG",
  "ATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGATCAGA"
  };

  for ( int j(1) ; j <= 5 ; ++j )
  {

    SequenceReaderString readerRepeats(repeats[j-1],cout);
    HashTable hashRepeats(cout);
    creator.createHashTable(hashRepeats,readerRepeats,10,20);

    readerRepeats.rewind();

    WordSequence seqRepeats;

    readerRepeats.getNextSequence( seqRepeats,10 );

    HitListVector hlFwd, hlRev;

    //    hlFwd.clear();
    //    hlRev.clear();

    cout << "hello\n";    
    
    hashRepeats.setNumRepeats(j);
    //    hashRepeats.matchSequence( seqRepeats, hlFwd, j );  
    hashRepeats.matchSequence( seqRepeats, hlFwd );  

    WordSequence rc;

    reverseComplement( seqRepeats, rc, 10 );
    
    //    hashRepeats.matchSequence( rc, hlRev, j );  
    hashRepeats.matchSequence( rc, hlRev );  


    cout << "Repeats masked: " << j 
         << "\nHits found in forward direction: " << hlFwd.size() 
         << "\nHits found in reverse direction: " << hlRev.size() << "\n";

    assert ( hlFwd.size() == 9 );

    if ( ( j % 2 ) == 0 ) assert( hlRev.size() == 9 );
    else assert( hlRev.size() == 0 );


      //    for ( HitListVector::iterator i( hlFwd.begin() ) ; i != hlFwd.end() ; ++i )
      //   cout << i->subjectNum << " " << i->diff << " " << i->queryPos << endl;

      //    cout << "rc\n";
      //   for ( HitListVector::iterator i( hlRev.begin() ) ; i != hlRev.end() ; ++i )
      //   cout << i->subjectNum << " " << i->diff << " " << i->queryPos << endl;


 
  } // ~for j

  for ( int j(1) ; j <= 5 ; ++j )
  {

    SequenceReaderString readerRepeats(repeats[j-1],cout);
    HashTablePacked hashRepeats(cout);
    creator.createHashTable(hashRepeats,readerRepeats,10,20);

    readerRepeats.rewind();

    WordSequence seqRepeats;

    readerRepeats.getNextSequence( seqRepeats,10 );

    HitListVector hlFwd, hlRev;

    //    hlFwd.clear();
    //    hlRev.clear();

    cout << "hello\n";    
    
    hashRepeats.setNumRepeats(j);
    //    hashRepeats.matchSequence( seqRepeats, hlFwd, j );  
    hashRepeats.matchSequence( seqRepeats, hlFwd );  

    WordSequence rc;

    reverseComplement( seqRepeats, rc, 10 );
    
    //    hashRepeats.matchSequence( rc, hlRev, j );  
    hashRepeats.matchSequence( rc, hlRev );  


    cout << "Repeats masked: " << j 
         << "\nHits found in forward direction: " << hlFwd.size() 
         << "\nHits found in reverse direction: " << hlRev.size() << "\n";

    assert ( hlFwd.size() == 9 );

    if ( ( j % 2 ) == 0 ) assert( hlRev.size() == 9 );
    else assert( hlRev.size() == 0 );


  } // ~for j

  
  cout << "Test passed!" << endl << endl; 

  // ---

  cout << "Test " << ++numTests 
       <<": test of HashTableTranslated\n\n"; 
  {
    SequenceReaderFasta testReader( "test.fasta" );
   
    int wordLength(4);

    HashTableTranslated hashTrans(cout,"test_trans");

    hashTrans.createHashTable( testReader, wordLength, 10000 );
    assert(hashTrans.isInitialized());

    hashTrans.saveHashTable();

    HashTableTranslated hashTrans2(cout,"test_trans");
    hashTrans2.loadHashTable();
    assert( hashTrans2.isInitialized() );
    assert( hashTrans.getNumSequences() == hashTrans2.getNumSequences() );
    string s1,s2;

    cout << hashTrans.getMaxNumHits() << " HH " << hashTrans2.getMaxNumHits() << endl;


    testReader.rewind();
    HitListVector h1, h2;
    WordSequence w1,w2;
    WordSequence t1,t2; // to store translations of w1, w2
    CodonList codons;
    SequenceEncoderCodon encoder;
    encoder.setWordLength(wordLength);

    for ( int i(1) ; i < hashTrans.getNumSequences() ; i++ )
    {

      hashTrans.getSequenceName( s1, i );
      hashTrans2.getSequenceName( s2, i );
      assert (s1 == s2 );
      assert (s1 == (string) hashTrans.getSequenceName(i) );
      assert (s1 == (string) hashTrans2.getSequenceName(i) );

      h1.clear(); h2.clear(); w1.clear(); w2.clear();
      testReader.getNextSequence( w1, gMaxBasesPerWord );
      cout << printWord( w1, wordLength ) << endl;
      w2=w1;

      assert(hashTrans.getSequenceSize(i)==hashTrans2.getSequenceSize(i));

      assert( (((w1.size()-1) * gMaxBasesPerWord ) + w1.getNumBasesInLast())
              == hashTrans.getSequenceSize(i));


      // Sequence should produce same (nonzero) num hits in fwd direction

      hashTrans.setForward(); assert(hashTrans.isForward());
      hashTrans2.setForward(); assert(hashTrans2.isForward());

      hashTrans.setQueryTranslatedDNA();
      hashTrans2.setQueryTranslatedDNA();
      hashTrans.matchSequence( w1, h1 );
      hashTrans2.matchSequence( w2, h2 );
      assert(h1.size()!=0);
      assert(h1.size()==h2.size());
      assert(h1==h2);
      assert(w1==w2); // translatedDNA matching does not affect w1 or w2

      // now encode to codons for each frame
      hashTrans.setQueryProtein();
      hashTrans2.setQueryProtein();
      for ( int rf(0) ; rf < gNumReadingFrames ; rf++ )
      {
        codons.clear();
        codonize( w1, codons, rf );
        encoder.linkSeq(t1);
        encoder.encode( codons );
        encoder.unlinkSeq();
        //      cout << printResidue( t1, wordLength ) << endl;
        t2=t1;
        h1.clear();h2.clear();
        hashTrans.matchSequence( t1, h1 );
        hashTrans2.matchSequence( t2, h2 );
        assert(h1.size()!=0);
        assert(h1.size()==h2.size());
        assert(h1==h2);
        
      } // ~for rf


      // RC of sequence should produce same (nonzero) num hits in rev direction

      w1.clear();
      reverseComplement( w2,w1, gMaxBasesPerWord );
      w2=w1;

      h1.clear(); h2.clear();
      hashTrans.setQueryTranslatedDNA();
      hashTrans2.setQueryTranslatedDNA();
      hashTrans.setReverse(); assert(!hashTrans.isForward());
      hashTrans2.setReverse(); assert(!hashTrans2.isForward());
      hashTrans.matchSequenceTranslatedDNA( w1, h1 );
      hashTrans2.matchSequenceTranslatedDNA( w2, h2 );
      assert(h1.size()!=0);
      assert(h1.size()==h2.size());
      assert(h1==h2);
      assert(w1==w2); // translatedDNA matching does not affect w1 or w2

      // now encode to codons for each frame
      hashTrans.setQueryProtein();
      hashTrans2.setQueryProtein();
      for ( int rf(0) ; rf < gNumReadingFrames ; rf++ )
      {
        codons.clear();
        codonize( w1, codons, rf );
        encoder.linkSeq(t1);
        encoder.encode( codons );
        encoder.unlinkSeq();
        //      cout << printResidue( t1, wordLength ) << endl;
        t2=t1;
        h1.clear();h2.clear();
        hashTrans.matchSequence( t1, h1 );
        hashTrans2.matchSequence( t2, h2 );
        assert(h1.size()!=0);
        assert(h1.size()==h2.size());
        assert(h1==h2);
        
      } // ~for rf

     
    } //~ for i


  } // ~test

  cout << "Test passed!\n\n"; 


  // ---

  cout << "Test " << ++numTests 
       <<": test of substitute generation\n\n";

  {
    string s, s1;
    Word w;
    vector<Word> subs;
    int wl=15;

    // test substitution for DNA
    // 

    for (int i(0); i<wl; i++) s+="A";

    w = makeBase(s);

    generateSubstitutesDNA( w, subs, wl );
    assert(subs.size()==wl);
    
    for (int i(0); i<wl ; i++)
    {
      cout << printWord(subs[i],wl) << endl;
      s1=s;
      s1[wl-i-1]='G';
      assert(makeBase(s1)==subs[i]);
    }

    //
    s=""; subs.clear();

    for (int i(0); i<wl; i++) s+="C";

    w = makeBase(s);

    generateSubstitutesDNA( w, subs, wl );
    assert(subs.size()==wl);
    
    for (int i(0); i<wl ; i++)
    {
      cout << printWord(subs[i],wl) << endl;
      s1=s;
      s1[wl-i-1]='T';
      assert(makeBase(s1)==subs[i]);
    }

    //
    s=""; subs.clear();

    for (int i(0); i<wl; i++) s+="G";

    w = makeBase(s);

    generateSubstitutesDNA( w, subs, wl );
    assert(subs.size()==wl);
    
    for (int i(0); i<wl ; i++)
    {
      cout << printWord(subs[i],wl) << endl;
      s1=s;
      s1[wl-i-1]='A';
      assert(makeBase(s1)==subs[i]);
    }

    //
    s=""; subs.clear();

    for (int i(0); i<wl; i++) s+="T";

    w = makeBase(s);

    generateSubstitutesDNA( w, subs, wl );
    assert(subs.size()==wl);
    
    for (int i(0); i<wl ; i++)
    {
      cout << printWord(subs[i],wl) << endl;
      s1=s;
      s1[wl-i-1]='C';
      assert(makeBase(s1)==subs[i]);
    }

    // test substitution values for protein
    for (int i(0); i<gNumCodonEncodings; i++)
    {
      s=""; s+= gResidueNames[i];
      w = makeResidue(s);
      subs.clear();
      generateSubstitutesProtein( w, subs, 1 );
      cout << "Residue: " << s << " substitutes: ";
        for (vector<Word>::iterator j(subs.begin());j!=subs.end();j++)
          cout << printResidue( *j, 1 ) << " ";
      cout << endl;
    }

    s="*CPXY";
    w=makeResidue(s);
    subs.clear();
    generateSubstitutesProtein( w, subs, 5 );
    for (vector<Word>::iterator j(subs.begin());j!=subs.end();j++)
      cout << printResidue( *j, 5 ) << endl;
    assert(subs.size()==3);

    s="EIKLMN";
    w=makeResidue(s);
    subs.clear();
    generateSubstitutesProtein( w, subs, 6 );
    for (vector<Word>::iterator j(subs.begin());j!=subs.end();j++)
      cout << printResidue( *j, 6 ) << endl;
    assert(subs.size()==18);


    assert(subVals[subStarts[gNumCodonEncodings]]==99999);

    // Test substitution for a DNA hash table
    {
    string orig = 
    "ttcgtgacaggcagttcaggtcacgtggtgacttgatgacccattgtcaaatgtccaattgtta";
    string subd =
    "CtcgtgacaAgcagttcaAgtcacgtgAtgacttgaCgacccattAtcaaatgtTcaattgttG";
 //  1234567812345678123456781234567812345678123456781234567812345678
    
    SequenceReaderString origReader(orig,cout);
    SequenceReaderString subdReader(subd,cout);

    HashTablePacked   origHash(cout);
    creator.createHashTable(origHash,origReader,8,100000);

    WordSequence subdSeq, seq;

    subdReader.getNextSequence( subdSeq,8 );

    HitListVector hl;

    seq=subdSeq;
    origHash.matchSequence( seq, hl );  
    assert(hl.size()==0);

    origHash.setSubstituteThreshold(1);
    
    seq=subdSeq;
    origHash.matchSequence( seq, hl );  
    for (HitListVector::iterator i(hl.begin()); i!=hl.end();i++)
      cout << i->subjectNum << " " << i->diff << " " << i->queryPos << endl;
    assert(hl.size()==8);
    }

    {
      string orig = 
        "LAKAGKKTGPQDAECRSCHNPKPSAASSWKEIGFDKSLHYRHVASKAIKPVGDPQKNCGA";
      string subd =
        "iAKAGKeTGPQDsECRSCnNPKPSsASSfKEIGFnKtLHYkHVASKAmKPVGDPkKhCGA";
//       123451234512345123451234512345123451234512345123451234512345
          
    SequenceReaderStringProtein origReader(orig,cout);
    SequenceReaderStringProtein subdReader(subd,cout);

    HashTablePackedProtein origHash(cout);
    creator.createHashTable(origHash,origReader,5,100000);

    WordSequence subdSeq, seq;

    subdReader.getNextSequence( subdSeq,5 );

    HitListVector hl;

    seq=subdSeq;
    origHash.matchSequence( seq, hl );  
    assert(hl.size()==0);

    origHash.setSubstituteThreshold(1);
    
    seq=subdSeq;
    origHash.matchSequence( seq, hl );  
    for (HitListVector::iterator i(hl.begin()); i!=hl.end();i++)
      cout << i->subjectNum << " " << i->diff << " " << i->queryPos << endl;
    cout << hl.size() << endl;
    assert(hl.size()==12);
    }


  } // ~test 
  cout << "Test passed!\n";

  // ---

  cout << "Test " << ++numTests 
       <<": test of flag mode\n\n"; 
  {
    string testSub
    ="TTTTTTT-TTTTTTT-TTTTTTT-TTTTTTT-TTTTTTT-TTTTTTT-TTTTTTT-TTTTTTT";
    //123456712345671234567123456712345671234567123456712345671234567
    string testQuery="TTTNTTT";
   
    HitListVector hl1, hl2;
    
    SequenceReaderString 
      queryReader(testQuery, cout), subjectReader(testSub,cout);
    SequenceReaderModeReplace replace('T');      
    SequenceReaderModeFlagReplace tag('T');
    WordSequence testSeq, seq;
    queryReader.changeMode(&replace);
    queryReader.getNextSequence( testSeq, 7 );
  



    {    
      subjectReader.changeMode(&replace);
      
      HashTable hash(cout);
      HashTablePacked hashPacked(cout);
      creator.createHashTable(hash,subjectReader,7,1000);
      creator.createHashTable(hashPacked,subjectReader,7,1000);

      seq=testSeq;
      hash.matchSequence( seq, hl1 );  
      for (HitListVector::iterator i(hl1.begin()); i!=hl1.end();i++)
        cout << i->subjectNum 
             << " " << i->diff << " " << i->queryPos << endl;
      cout << hl1.size() << endl;
      assert(hl1.size()==9);

      seq=testSeq;
      hashPacked.matchSequence( seq, hl2 );  
      assert(hl1==hl2);

      hl1.clear();
      hl2.clear();

      // Tagging of subject not implemented properly %%%%%%%
      // check tagging of query works
      queryReader.changeMode(&tag);
      queryReader.rewind();
      queryReader.getNextSequence( seq, 7 );
      hash.matchSequence( seq, hl1 );  
      hashPacked.matchSequence( seq, hl2 );  
      assert(hl1.size()==0);
      assert(hl2.size()==0);

    }
    {
      // now check tagging of subject works
      subjectReader.changeMode(&tag);
      hl1.clear();hl2.clear();
      
      HashTable hash(cout);
      HashTablePacked hashPacked(cout);
      creator.createHashTable(hash,subjectReader,7,1000);
      creator.createHashTable(hashPacked,subjectReader,7,1000);

      seq=testSeq;
      hash.matchSequence( seq, hl1 );  

      for (HitListVector::iterator i(hl1.begin()); i!=hl1.end();i++)
        cout << i->subjectNum 
             << " " << i->diff << " " << i->queryPos << endl;
      cout << hl1.size() << endl;
      assert(hl1.size()==2);

      seq=testSeq;
      hashPacked.matchSequence( seq, hl2 );  

      for (HitListVector::iterator i(hl2.begin()); i!=hl2.end();i++)
        cout << i->subjectNum 
             << " " << i->diff << " " << i->queryPos << endl;
      cout << hl2.size() << endl;

      assert(hl1==hl2);


    }


  } // ~test
  cout << "Test passed!" << endl << endl; 


  // ---

  cout << "Test " << ++numTests 
       <<": test of flag mode for translated DNA\n\n"; 
  {
    
    SequenceReaderModeReplace replaceMode('A');
    SequenceReaderModeFlagReplace flagMode('A');
    string testSeq = 
       "agtagtagtagtagtagtagtagtagtagtXgtagtagtagtagtagtagtagtagtagtag";
    //  1..|..|..|..|..2..|..|..|..|..3..|..|..|..|..4..|..|..|..|..
    // Forward direction: f0=agt=SSS..., f1=gta=VVV....., f2=tag=***...
    // Reverse direction: r0=act=TTT..., r1=cta=LLL....., r2=tac=YYY...
    SequenceReaderString subjectReader(testSeq, cout);

    HashTableTranslatedPeek tableUnflagged, tableFlagged;

    HashTableComponent& peekUnflaggedFwd(tableUnflagged.peekFwd());
    HashTableComponent& peekUnflaggedRev(tableUnflagged.peekRev());
    HashTableComponent& peekFlaggedFwd(tableFlagged.peekFwd());
    HashTableComponent& peekFlaggedRev(tableFlagged.peekRev());

    cout << "making unflagged" << endl;
    subjectReader.changeMode(&replaceMode);
    tableUnflagged.createHashTable( subjectReader, 5, 10000 );
    cout << endl << "  **** making flagged ****" << endl << endl;
    subjectReader.changeMode(&flagMode);
    tableFlagged.createHashTable( subjectReader, 5, 10000 );

    string queryString="SSSSSSSSS";

    SequenceReaderStringProtein queryReader(queryString,cout);
    WordSequence w1, w2;
    HitListVector h1, h2;
    queryReader.getNextSequence(w1,5);
    assert(w1.size()==2);
    //    assert(w1.getNumBasesInLast()==0);
    w2=w1;

    cout << "hits for unflagged" << endl;

    tableUnflagged.matchSequenceProtein( w1, h1);
    for (HitListVector::iterator i(h1.begin()); i!=h1.end(); i++)
      cout << i->subjectNum << " " << i->diff << " " << i->queryPos << endl;

    cout << "hits for flagged" << endl;

    tableFlagged.matchSequenceProtein( w2, h2);
    for (HitListVector::iterator i(h2.begin()); i!=h2.end(); i++)
      cout << i->subjectNum << " " << i->diff << " " << i->queryPos << endl;




    Word w = makeResidue((string)"SSSSS");
    PackedHitStore phu, phf;
    HitListVector hu, hf;

    cout << "matching unflagged" << endl;
    peekUnflaggedFwd.matchWordStandard(w, phu, 0);
    peekUnflaggedFwd.convertHits( phu, hu );

    for (HitListVector::iterator i(hu.begin()); i!=hu.end(); i++)
      cout << i->subjectNum << " " << i->diff << " " << i->queryPos<< endl;

    cout << "matching flagged" << endl;
    peekFlaggedFwd.matchWord(w, phf, 0);
    cout << "converting hits" << endl;
    peekFlaggedFwd.convertHits( phf, hf );

    for (HitListVector::iterator i(hf.begin()); i!=hf.end(); i++)
      cout << i->subjectNum << " " << i->diff << " " << i->queryPos << endl;

    hu.clear(); hf.clear(); phu.clear(); phf.clear();
    
    w=makeResidue((string)"YYYYY");


    cout << "matching unflagged" << endl;
    peekUnflaggedRev.matchWordStandard(w, phu, 0);
    peekUnflaggedRev.convertHits( phu, hu );

    for (HitListVector::iterator i(hu.begin()); i!=hu.end(); i++)
      cout << i->subjectNum << " " << i->diff << " " << i->queryPos<< endl;

    cout << "matching flagged" << endl;
    peekFlaggedRev.matchWord(w, phf, 0);
    cout << "converting hits" << endl;
    peekFlaggedRev.convertHits( phf, hf );

    for (HitListVector::iterator i(hf.begin()); i!=hf.end(); i++)
      cout << i->subjectNum << " " << i->diff << " " << i->queryPos << endl;


  } // ~test 
  cout << "Test passed!\n";



  return (0);

  } // ~try
  catch (const SSAHAException& err )
  {
    cout << "Caught SSAHA exception: " << err.what() << "\n";
    exit(1);
  }
  catch (const std::exception& err )
  {
    cout << "Caught exception: " << err.what() << "\n";
    exit(1);
  }


} // ~main for testHashTable

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

// End of file testHashTable.cpp

---