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

// Description:

// Includes:

#include <exception>
#include "SequenceEncoder.h"
#include "SequenceReader.h"
#include "SequenceReaderFasta.h"
#include "SequenceReaderFastq.h"
#include "SequenceReaderString.h"
#include "SequenceReaderMulti.h"
#include "SequenceReaderLocal.h"
#include "SequenceReaderFilter.h"
//#include "SequenceReaderCodon.h"
#include "GenerateTestFastaFiles.h"
#include "assert.h"
#include <string>
#include <strstream>
#include <iostream>

// ### Function Definitions ###
void capitalise( string& s )
{
  for ( int i(0) ; i < s.size() ; ++i ) s[i] = toupper(s[i]);
}
void reverseString( string& seq )
{
  string rc;
  for ( int i(0) ; i < seq.size() ; i++ )
  {
    if  ( ( seq[i] == 'A' ) || (seq[i] == 'a') ) rc = 'T' + rc; 
    else if ( ( seq[i] == 'T' ) || (seq[i] == 't') ) rc = 'A' + rc; 
    else if ( ( seq[i] == 'G' ) || (seq[i] == 'g') ) rc = 'C' + rc; 
    else if ( ( seq[i] == 'C' ) || (seq[i] == 'c') ) rc = 'G' + rc; 
    else throw SSAHAException();
  } // ~for i

  seq = rc;

} // ~reverseString




int main( void )
{
  try
  {

  int numTests = 0;

  cout << "*******************************************" << endl << endl;
  cout << "Test of functions in SequenceReader modules" << endl << endl;
  cout << "*******************************************" << endl << endl;

  const int wordLength(10);
  std::ostrstream buffer;
  string s1,s2,s3,s4,s5;

  // ---

  cout << " --- Tests of member functions of SequenceReaderFasta ---\n\n";
  cout << "Test " << ++numTests <<": test of function getNextSequence" 
       << endl << endl;

  int numSeqs = 10;
  int seqSize = 57;

  // Generate a random sequence of (numSeqs*seqSize) base pairs ...
  BaseGenerator testBases( numSeqs * seqSize );

  // ... 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);
  WordSequence w;
  string expected, actual, name;

  // last Word of decoded WordSequence will not be full, and the blank spaces
  // will be decoded by printWord as 'A's. spareBases is a string of 'A's to
  // add on the end of the expected sequence string for it to match actual.
  //  const string spareBases( wordLength - ( seqSize % wordLength ), 'A');

  // TC 08.03.01 - improved printWord so it takes account of the number of
  // bases in the last word of a WordSequence (by looking at numBasesInLast)
  // Therefore above step no longer necessary.
  const string spareBases("");



  for ( int i(0) ; i < numSeqs ; i++ )
  {
  
    testReader.getNextSequence(w,wordLength);
    buffer << printWord(w,wordLength) << ends;
    actual=buffer.str();
    buffer.freeze(false); // release memory to write further stuff to buffer  
    buffer.seekp(0,ios::beg); // ensure next write is to start of buffer
 
    testBases.getBases(i*seqSize,seqSize,expected);

    expected += spareBases;

    testReader.getLastSequenceName(name);
    cout << "SequenceName     : " << name << endl;
    cout << "Expected sequence: " << expected << endl;
    cout << "Actual sequence  : " << actual   << endl;

    assert(expected == actual);

    w.clear();

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

  // ---

  cout << "Test " << ++numTests <<": test of sequence numbering" 
       << endl << endl;


  cout << "Last sequence number read: " << testReader.getLastSequenceNumber()
       << endl;

  assert( testReader.getLastSequenceNumber() == numSeqs );

  int numInLast = testReader.getNextSequence( w, wordLength );

  cout << "getNextSequence returned " << numInLast 
       << ": should be -1 because all sequences have been read." << endl;

  assert( numInLast == -1 );

  cout << "Number of sequences in file: " << testReader.getNumSequencesInFile()
       << endl;

  assert( testReader.getNumSequencesInFile() == numSeqs );

  cout << "Rewinding to start of data ..." << endl;

  testReader.rewind();

  cout << "Last sequence number read: " << testReader.getLastSequenceNumber()
       << endl;

  assert( testReader.getLastSequenceNumber() == 0 );
 
  cout << "Checking getBitsPerSymbol returns 2" << endl;

  assert( testReader.getBitsPerSymbol()==2 );

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


  // ---

  cout << "Test " << ++numTests <<": test of getSequence" 
       << endl << endl;

    int toRead = 5;

    w.clear();
    testReader.getSequence(w,toRead,wordLength);
    buffer << printWord(w,wordLength) << ends;
    actual=buffer.str();
    buffer.freeze(false); // release memory to write further stuff to buffer  
    buffer.seekp(0,ios::beg); // ensure next write is to start of buffer
 
    testBases.getBases((toRead-1)*seqSize,seqSize,expected);

    expected += spareBases;

    cout << "Expected sequence: " << expected << endl;
    cout << "Actual sequence  : " << actual   << endl;

    assert(expected == actual);

    cout << toRead << " " << testReader.getLastSequenceNumber() << endl;
    assert(testReader.getLastSequenceNumber() == toRead );    

    w.clear();


    //    cout << 
    //   "Passing invalid sequence num to getSequence, should throw exception...\n";

    //    try
    //   {
    //    numInLast = testReader.getSequence( w, 987, wordLength );
    //    cout << "Error: exception should have been thrown!\n";
    //     exit(1);
    //   } // ~try
    //   catch (const SSAHAException& err )
    //   {
    //     cout << "Caught SSAHA exception: " << err.what() << "\n";
    //   }

    cout << "Passing invalid sequence number to getSequence ..." << endl;

    numInLast = testReader.getSequence( w, 987, wordLength );

    cout << "getNextSequence returned " << numInLast 
         << ": should be -1 because sequence number is invalid." << endl;

    assert( numInLast == -1 );

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

    // ---

    cout << "Test " << ++numTests <<": test of random access output functions" 
       << endl << endl;

    for ( int i(1) ; i <= numSeqs; i++ )
    {
      cout << testReader.getName(i) << endl;
      cout << testReader.getSideInfo(i) << endl;
      cout << testReader.getSource(i) << endl;
    } // ~for i

  // ---

  cout << " --- Tests of member functions of SequenceReaderFastq ---\n\n";
  cout << "Test " << ++numTests <<": test of function getNextSequence" 
       << endl << endl;

  //  int numSeqs = 10;
  // int seqSize = 57;

  // Generate a random sequence of (numSeqs*seqSize) base pairs ...
  //  BaseGenerator testBases( numSeqs * seqSize );

  // ... 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.generateSubjectFileFastq(seqSize,0);

  SequenceReaderFastq testReaderFastq("test_subject.fastq",cout);
  //  WordSequence w;
  //  string expected, actual, name;

  // last Word of decoded WordSequence will not be full, and the blank spaces
  // will be decoded by printWord as 'A's. spareBases is a string of 'A's to
  // add on the end of the expected sequence string for it to match actual.
  //  const string spareBases( wordLength - ( seqSize % wordLength ), 'A');

  for ( int i(0) ; i < numSeqs ; i++ )
  {
  
    testReaderFastq.getNextSequence(w,wordLength);
    buffer << printWord(w,wordLength) << ends;
    actual=buffer.str();
    buffer.freeze(false); // release memory to write further stuff to buffer  
    buffer.seekp(0,ios::beg); // ensure next write is to start of buffer
 
    testBases.getBases(i*seqSize,seqSize,expected);

    expected += spareBases;

    testReaderFastq.getLastSequenceName(name);
    cout << "SequenceName     : " << name << endl;
    cout << "Expected sequence: " << expected << endl;
    cout << "Actual sequence  : " << actual   << endl;

    assert(expected == actual);

    w.clear();

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

  // ---

  cout << "Test " << ++numTests <<": test of sequence numbering" 
       << endl << endl;


  cout << "Last sequence number read: " << testReaderFastq.getLastSequenceNumber()
       << endl;

  assert( testReaderFastq.getLastSequenceNumber() == numSeqs );

  numInLast = testReaderFastq.getNextSequence( w, wordLength );

  cout << "getNextSequence returned " << numInLast 
       << ": should be -1 because all sequences have been read." << endl;

  assert( numInLast == -1 );

  cout << "Number of sequences in file: " << testReaderFastq.getNumSequencesInFile()
       << endl;

  assert( testReaderFastq.getNumSequencesInFile() == numSeqs );

  cout << "Rewinding to start of data ..." << endl;

  testReaderFastq.rewind();

  cout << "Last sequence number read: " << testReaderFastq.getLastSequenceNumber()
       << endl;

  assert( testReaderFastq.getLastSequenceNumber() == 0 );
 
  cout << "Test passed!" << endl << endl; 

  cout << "Checking getBitsPerSymbol returns 2" << endl;

  assert( testReader.getBitsPerSymbol()==2 );

  // ---

  cout << "Test " << ++numTests <<": test of getSequence" 
       << endl << endl;

    toRead = 5;

    w.clear();
    testReaderFastq.getSequence(w,toRead,wordLength);
    buffer << printWord(w,wordLength) << ends;
    actual=buffer.str();
    buffer.freeze(false); // release memory to write further stuff to buffer  
    buffer.seekp(0,ios::beg); // ensure next write is to start of buffer
 
    testBases.getBases((toRead-1)*seqSize,seqSize,expected);

    expected += spareBases;

    cout << "Expected sequence: " << expected << endl;
    cout << "Actual sequence  : " << actual   << endl;


    assert(expected == actual);

    cout << toRead << " " << testReaderFastq.getLastSequenceNumber() << endl;
    assert(testReaderFastq.getLastSequenceNumber() == toRead );    

    w.clear();

    //    cout << 
    //   "Passing invalid sequence num to getSequence, should throw exception...\n";

    //    try
    //    {
    //    numInLast = testReaderFastq.getSequence( w, 987, wordLength );
    //     cout << "Error: exception should have been thrown!\n";
    //     exit(1);
    //   } // ~try
    //   catch (const SSAHAException& err )
    //   {
    //     cout << "Caught SSAHA exception: " << err.what() << "\n";
    //   }

    //    cout << "getNextSequence returned " << numInLast 
    //        << ": should be -1 because sequence number is invalid." << endl;

    //    assert( numInLast == -1 );

    cout << "Passing invalid sequence number to getSequence ..." << endl;

    numInLast = testReaderFastq.getSequence( w, 987, wordLength );

    cout << "getNextSequence returned " << numInLast 
         << ": should be -1 because sequence number is invalid." << endl;

    assert( numInLast == -1 );

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

    // ---

    cout << "Test " << ++numTests <<": test of random access output functions" 
       << endl << endl;

    for ( int i(1) ; i <= numSeqs; i++ )
    {
      cout << testReaderFastq.getName(i) << endl;
      cout << testReaderFastq.getSideInfo(i) << endl;
      cout << testReaderFastq.getSource(i) << endl;
    } // ~for i

    // ---

    cout << "Test " << ++numTests <<": test of SequenceReaderString" 
         << endl << endl;
    //        12345678901234567890123456789012345678901234567890
    expected="ATGCCTGCGTATGTAATGCTGGCATTGGTGCGTTTAATGCTCTCCTGATG";
    SequenceReaderString stringReader(expected,cout);
    expected+="AAAAAAAAAA"; //%%%

    assert(stringReader.getNumSequencesInFile()==1);

    w.clear();
    numInLast = stringReader.getNextSequence(w,wordLength);
    buffer << printWord(w,wordLength) << ends;
    actual=buffer.str();
    buffer.freeze(false); // release memory to write further stuff to buffer  
    buffer.seekp(0,ios::beg); // ensure next write is to start of buffer
    
    cout << "expected: " << expected << endl;
    cout << "actual:   " << actual << " " << w.size() << " " 
         << w.getNumBasesInLast() << endl;


    assert( numInLast == 0 );
    assert( expected == actual );    

    cout << "Checking getBitsPerSymbol returns 2" << endl;

    assert( stringReader.getBitsPerSymbol()==2 );



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

    // ---

    cout << "Test " << ++numTests <<": test of SequenceReaderMulti" 
         << endl << endl;

    SequenceReaderMulti multiReader(cout);

    int numReaders = 4;
    
    for ( int i(0) ; i < numReaders ; i ++ ) 
      multiReader.addReader(testReader.clone());

    assert( multiReader.getNumSequencesInFile() == 
            ( numReaders * testReader.getNumSequencesInFile() ) );

    WordSequence wSingle, wMulti;
    multiReader.rewind();

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

      testReader.rewind();
      for ( int j(0) ; j < testReader.getNumSequencesInFile() ; j ++ ) 
      {
                cout << i << " " << j << endl;
        wSingle.clear(); wMulti.clear(); 
        assert( multiReader.getBitsPerSymbol()==2 );
        assert( testReader.getNextSequence(wSingle, wordLength) ==  
                multiReader.getNextSequence(wMulti, wordLength) );
        //         cout << multiReader.getLastSequenceNumber() << " " 
        //              << testReader.getLastSequenceNumber() << endl;
        testReader.getLastSequenceName(expected);
        multiReader.getLastSequenceName(actual);
        cout << expected << " " << actual << endl;
        assert( expected == actual );

        assert( multiReader.getLastSequenceNumber()
                == ( i*testReader.getNumSequencesInFile() )
                   + testReader.getLastSequenceNumber() );
        assert ( wSingle == wMulti );
      } // ~for j

    } // ~for i

    wMulti.clear();
   
    cout << multiReader.getLastSequenceNumber()  << endl;

    assert( multiReader.getLastSequenceNumber() == 
            ( numReaders * testReader.getNumSequencesInFile() ) );
    assert( multiReader.getNextSequence( wMulti,wordLength ) == -1 );

    multiReader.rewind();

    assert( multiReader.getLastSequenceNumber() == 0 );   

    cout << "Checking getBitsPerSymbol returns 2" << endl;

    assert( multiReader.getBitsPerSymbol()==2 );

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

    // ---

    cout << "Test " << ++numTests <<": test of make/printBase/Residue" 
         << endl << endl;

    {
      Word w1;
      string s1;

      s1="TaCgTTgtAGc"; 
      w1=makeBase(s1);
      buffer << printBase(w1,11) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << actual << "\n";
      cout << s1 << "\n";
      capitalise(s1);
      assert( actual == s1 );

      s1="ACDEFG"; 
      w1=makeResidue(s1);
      buffer << printResidue(w1,6) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << actual << "\n";
      cout << s1 << "\n";
      capitalise(s1);
      assert( actual == s1 );

      s1="NMLKIH"; 
      w1=makeResidue(s1);
      buffer << printResidue(w1,6) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << actual << "\n";
      cout << s1 << "\n";
      capitalise(s1);
      assert( actual == s1 );

      //      s1="PRTQSU"; 
      s1="PRTQSX"; 
      w1=makeResidue(s1);
      buffer << printResidue(w1,6) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << actual << "\n";
      cout << s1 << "\n";
      capitalise(s1);
      assert( actual == s1 );

      s1="VWY***"; 
      w1=makeResidue(s1);
      buffer << printResidue(w1,6) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << actual << "\n";
      cout << s1 << "\n";
      capitalise(s1);
      assert( actual == s1 );

      s1="acdefg"; 
      w1=makeResidue(s1);
      buffer << printResidue(w1,6) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << actual << "\n";
      cout << s1 << "\n";
      capitalise(s1);
      assert( actual == s1 );

      s1="nmlkih"; 
      w1=makeResidue(s1);
      buffer << printResidue(w1,6) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << actual << "\n";
      cout << s1 << "\n";
      capitalise(s1);
      assert( actual == s1 );

      //      s1="prtqsu"; 
      s1="prtqsx"; 
      w1=makeResidue(s1);
      buffer << printResidue(w1,6) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << actual << "\n";
      cout << s1 << "\n";
      capitalise(s1);
      assert( actual == s1 );

      s1="vwy***"; 
      w1=makeResidue(s1);
      buffer << printResidue(w1,6) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << actual << "\n";
      cout << s1 << "\n";
      capitalise(s1);
      assert( actual == s1 );


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


    // ---

    cout << "Test " << ++numTests <<": test of SequenceEncoder"
         << endl << endl;
    {
      WordSequence W;
      SequenceEncoderDNA encoder(12);

      cout << "Checking getBitsPerSymbol returns 2" << endl;

      assert( encoder.getBitsPerSymbol()==2 );

      cout << "Checking mode changes...\n";

      string s1="nOne@rEv@lID";

      cout << "... by default, should ignore invalid characters\n";
      encoder.linkSeq(W);
      encoder.encode(s1);
      encoder.unlinkSeq();
      cout << W.size() << " " << W[0] << endl;
      // puts a single word into seq, but does nothing to it
      assert(W.size()==1);
      assert(W[0]==0);

      cout << "...changing mode to replace invalid chars with 'G'\n";
      SequenceReaderModeReportReplace mode1('G');
      //      encoder.changeMode( new SequenceReaderModeReportReplace('G') );
      encoder.changeMode( &mode1 );
      encoder.linkSeq(W);
      encoder.encode(s1);
      encoder.unlinkSeq();
      cout << W.size() << " " << printBase(W[0],12) << endl;
      assert(W.size()==2);
      assert(W[0]==makeBase("GGGGGGGGGGGG"));

      cout << "...checking mode changes are passed when copied\n";
      W.clear();
      SequenceEncoder copy(encoder);
      // shouldn't have to link, should be linked to same as encoder
      copy.linkSeq(W);
      copy.encode(s1);
      copy.unlinkSeq();
      cout << W.size() << " " << printBase(W[0],12) << endl;
      assert(W.size()==2);
      assert(W[0]==makeBase("GGGGGGGGGGGG"));

      cout << "...changing copy to 'ignore' mode\n";
      W.clear();
      SequenceReaderModeIgnore mode2;
      //      copy.changeMode( new SequenceReaderModeIgnore );
      copy.changeMode( &mode2 );
      copy.linkSeq(W);
      copy.encode(s1);
      copy.unlinkSeq();
      // puts a single word into seq, but does nothing to it
      assert(W.size()==1);
      assert(W[0]==0);

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


    // ---

    cout << "Test " << ++numTests <<": test of SequenceEncoderDNA" 
         << endl << endl;

    {
      //  cout << "for ttDNA\n";
      //  print(ttDNA);
      //  cout << "for ttProtein\n";
      //  print(ttProtein);
      //  cout << "for ttCodon\n";
      //  print(ttCodon);

      int wordLength(12);

      SequenceEncoderDNA encoder(wordLength);

      string s1="AGCcGTGGGTTTG", s2="TTGgGTTGGGGAC", s3="AGCTGaTTTGCAAGtCAA";

      WordSequence W;

      encoder.linkSeq(W);
      encoder.encode(s1);
      encoder.encode(s2);
      encoder.encode(s3);
      encoder.unlinkSeq();
      buffer << printBase(W,wordLength) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << printBase(W,wordLength) << "\n";
      cout << s1+s2+s3 << "\n";

      expected = s1+s2+s3;
      capitalise(expected);
      cout << actual << actual.size() << endl << expected << expected.size() << endl;

      assert(actual==expected);

      W.clear();

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

    // ---

    cout << "Test " << ++numTests <<": additional tests of SequenceEncoderDNA" 
         << endl << endl;

    {

      string fiftyAs="AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA";
      string thisString;
      WordSequence w;
      int numInLast;
      SequenceReaderModeIgnore ignore;
      SequenceReaderModeReplace replace('T');      
      SequenceReaderModeFlagReplace tag('T');
      int dudPos; 
      Word dudValue;

      
      for (int wordLength(9); wordLength <= 12; wordLength++)
      {
        for (int i(0) ; i < fiftyAs.size(); i++)
        {
          thisString=fiftyAs;
          thisString[i]='X';
          w.clear();
          dudPos=i/wordLength;
          dudValue=3<<(2*(wordLength-(i%wordLength)-1));
          cout << "dud " << wordLength << " " << i << " " 
               << dudPos << " " << printWord(dudValue, wordLength) << endl;

          SequenceReaderString stringReader(thisString,cout);


          assert(stringReader.getNumSequencesInFile()==1);

          // ignore mode
          cout << "testing ignore mode" << endl;

          stringReader.changeMode(&ignore);
          w.clear();
          numInLast = stringReader.getNextSequence(w,wordLength);
          assert(numInLast==(fiftyAs.size()-1)%wordLength);
          assert((w.size()-1)*wordLength+numInLast==fiftyAs.size()-1);

          for (WordSequence::iterator j(w.begin()); j!=w.end();j++)
          {
            assert (*j==0);
          }

          // replace mode
          cout << "testing replace mode" << endl;

          stringReader.rewind();
          stringReader.changeMode(&replace);
          w.clear();
          numInLast = stringReader.getNextSequence(w,wordLength);
          cout << "blegh " << wordLength << " " << i  << " " << fiftyAs.size() << " " << numInLast << " " << w.size() << endl;
          
          assert(numInLast==fiftyAs.size()%wordLength);
          assert((w.size()-1)*wordLength+numInLast==fiftyAs.size());

          for (WordSequence::iterator j(w.begin()); j!=w.end();j++)
          {
            cout << wordLength << " " << i << " " << (j-w.begin()) << " "
                 << *j << " " << gCursedWord << endl;
            if (j-w.begin()==dudPos) assert (*j==dudValue);
            else assert (*j==0);
          }

          // flag mode
          cout << "testing flag mode" << endl;

          stringReader.rewind();
          stringReader.changeMode(&tag);
          w.clear();
          numInLast = stringReader.getNextSequence(w,wordLength);
          assert(numInLast==fiftyAs.size()%wordLength);
          assert((w.size()-1)*wordLength+numInLast==fiftyAs.size());

          for (WordSequence::iterator j(w.begin()); j!=w.end();j++)
          {
            cout << wordLength << " " << i << " " << (j-w.begin()) << " "
                 << *j << " " << gCursedWord << endl;
            if (j-w.begin()==dudPos) assert (*j==dudValue|gCursedWord);
            else assert (*j==0);
          }
          


        } // ~for i


      } // ~test

      //  cout << "for ttDNA\n";
      //  print(ttDNA);
      //  cout << "for ttProtein\n";
      //  print(ttProtein);
      //  cout << "for ttCodon\n";
      //  print(ttCodon);

      int wordLength(12);

      SequenceEncoderDNA encoder(wordLength);

      string s1="AGCcGTGGGTTTG", s2="TTGgGTTGGGGAC", s3="AGCTGaTTTGCAAGtCAA";

      WordSequence W;

      encoder.linkSeq(W);
      encoder.encode(s1);
      encoder.encode(s2);
      encoder.encode(s3);
      encoder.unlinkSeq();
      buffer << printBase(W,wordLength) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << printBase(W,wordLength) << "\n";
      cout << s1+s2+s3 << "\n";

      expected = s1+s2+s3;
      capitalise(expected);
      cout << actual << actual.size() << endl << expected << expected.size() << endl;

      assert(actual==expected);

      W.clear();

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

    // ---

    cout << "Test " << ++numTests <<": test of SequenceReaderStringProtein" 
         << endl << endl;

    { //                 1....2....3....4....5....6....7....
      //      string protString("THEQUICK*R*WNF***UMPS*VERTHELA*YD*G");
      string protString("THEQXICK*R*WNF***XMPS*VERTHELA*YD*G");

      SequenceReaderStringProtein sr(protString,cout);

      int wordLength(5);
      WordSequence W, W2;
      sr.getNextSequence(W,wordLength);

      assert( sr.getNextSequence(W,wordLength) == -1 );

      sr.rewind();

      sr.getNextSequence(W2,wordLength);

      assert( W==W2 );

      string wordString;

      cout << "expected: " << protString << endl;
      cout << "actual:   " << printResidue(W,wordLength) << endl;
      

      for ( int i(0),j(0);  i < protString.size() ; i+=wordLength,j++ )
        {
          wordString = protString.substr(i, wordLength);
          cout << wordString << " " << printResidue(W[j],wordLength) << endl;
          assert(W[j]==makeResidue(wordString));
        }


      cout << "Checking getBitsPerSymbol returns 5" << endl;

      assert( sr.getBitsPerSymbol()==5 );

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

    // ---

    cout << "Test " << ++numTests <<": test of SequenceEncoderProtein" 
         << endl << endl;

    {

      int wordLength(5);

      SequenceEncoderProtein encoder(wordLength);

      cout << "Checking getBitsPerSymbol returns 5" << endl;

      assert( encoder.getBitsPerSymbol()==5 );

      //      string s1="YNLHEPCSADKFWQMV", s2="GI*TRUnvqramyg", s3="pswhcdfluteik";
      string s1="YNLHEPCSADKFWQMV", s2="GI*TRXnvqramyg", s3="pswhcdflxteik";

      WordSequence W;

      encoder.linkSeq(W);
      encoder.encode(s1);
      encoder.encode(s2);
      encoder.encode(s3);
      encoder.unlinkSeq();
      buffer << printResidue(W,wordLength) << ends;
      actual=buffer.str();
      buffer.freeze(false); 
      // release memory to write further stuff to buffer  
      buffer.seekp(0,ios::beg); 
      // ensure next write is to start of buffer
      cout << printResidue(W,wordLength) << "\n";
      cout << s1+s2+s3 << "\n";

      expected = s1+s2+s3; 
      capitalise(expected);
      assert(actual==expected);

      W.clear();

    }

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

    // ---

    cout << "Test " << ++numTests <<": test of codonize functions" 
         << endl << endl;

    {
      string allCodonsDNA ="taagcatgtgatgaatttggtcacataaaacttatgaatccccagcgttcaaccgtttggtat";

      CodonList c1, c2;
      WordSequence w1,w2,w1Rev;

      for (int k(0); k < 15; k++)
      {
        w1.clear(); w2.clear(); w1Rev.clear();
      SequenceReaderString reader(allCodonsDNA);
        reader.rewind();
        reader.getNextSequence(w1,gMaxBasesPerWord);
        reader.rewind();
        reader.getNextSequence(w2,gMaxBasesPerWord-1);
 
      for (int i(0) ; i <gNumReadingFrames; i++ )
      {
        c1.clear(); c2.clear();
        codonize( w1, c1, i );
        codonizeAndFlag( w2, c2, i );
        cout << "c1: " << c1 << endl;
        cout << "c2: " << c2 << endl;
        assert (c1==c2);
      }

      reverseComplement(w1, w1Rev, gMaxBasesPerWord);
      

      for (int i(0) ; i <gNumReadingFrames; i++ )
      {
        c1.clear(); c2.clear();
        codonize( w1Rev, c1, i );
        codonizeAndFlagReverse( w2, c2, i );
        cout << "c1: " << c1 << endl;
        cout << "c2: " << c2 << endl;
        assert (c1==c2);
      }

      allCodonsDNA=allCodonsDNA.substr(0,allCodonsDNA.size()-1);
      } // ~for k

    } // ~test

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


    // ---

#ifdef SEQUENCE_READER_CODON_NO_LONGER_USED
    cout << "Test " << ++numTests <<": test of SequenceReaderCodon" 
         << endl << endl;

    {


      int wordLength(5);
      const string expectedCodons((string)gCodonNames);
      SequenceReaderStringProtein exp(expectedCodons);
      WordSequence expectedSeq;
      exp.getNextSequence( expectedSeq, wordLength );
      cout << printResidue( expectedSeq, wordLength ) << endl;


      string allPossibleCodons;
      string allAs("");
      for ( int i(0) ; i < 4 ; i++ )
        for ( int j(0) ; j < 4 ; j ++ )
          for ( int k(0) ; k < 4 ; k ++ ) 
            {
              allPossibleCodons += gBaseNames[i];
              allPossibleCodons += gBaseNames[j];
              allPossibleCodons += gBaseNames[k];
            }


      cout << allPossibleCodons << endl;

      for ( int readFrame(0); readFrame < 6 ; readFrame ++ )
      {
        if (readFrame ==3) allAs="";
        string toBeCoded(allAs+allPossibleCodons);
        if (readFrame>=3) reverseString(toBeCoded);
     

        SequenceReaderString readDNA(toBeCoded,cout);
        SequenceReaderCodon  readCodon(readDNA,cout);

        assert( readCodon.getBitsPerSymbol()==5 );

        WordSequence w;
        string name;

        readDNA.getNextSequence( w, 12 );
        cout << printBase( w, 12 ) << endl;
        readDNA.rewind(); w.clear();

        for ( int perm(0) ; perm < 6 ; perm++ )
        {
          readCodon.getNextSequence( w, wordLength );
          readCodon.getLastSequenceName(name);
          cout << name << endl;
          cout << printResidue( w, wordLength ) << endl;        
          if (readFrame==perm) 
          {
            cout << printResidue( expectedSeq, wordLength ) << endl;    
            assert(w==expectedSeq);
          } // ~if
        } // ~for
        // shift read frame by 1
        allAs+="A";
        assert(readCodon.getNextSequence( w, wordLength )==-1);

      }

      cout << "testing SequenceEncoderCodon...\n";

 //       SequenceReaderString reader(allPossibleCodons);
//        WordSequence fd,rv, fdout, rvout;
//        reader.getNextSequence( fd, 16 );
//        reverseComplement( fd, rv, 16 );

//        SequenceEncoderCodon enc(5);
      
//        CodonList codons;

//        enc.setWordLength(5);

//        enc.linkSeq(fdout);

//        codonize( fd, codons, 0 );

//        enc.encode( codons);

//        enc.unlinkSeq();

//        codons.clear();

//        enc.setReverse();

//        enc.linkSeq(rvout);

//        codonize( rv, codons, 0 );

//        enc.encode( codons );
      
//        enc.unlinkSeq();

//        string sf,sr;

//        buffer << printResidue(fdout,5) << ends;
//        sf=buffer.str();
//        buffer.freeze(false); 
//        // release memory to write further stuff to buffer  
//        buffer.seekp(0,ios::beg); 

//        buffer << printResidue(rvout,5) << ends;
//        sr=buffer.str();
//        buffer.freeze(false); 
//        // release memory to write further stuff to buffer  
//        buffer.seekp(0,ios::beg); 

//        cout << sf << endl << sr << endl;

//        reverse( sr.begin(), sr.end() );

//        cout << sf << endl << sr << endl;

//        assert (sf==sr);



    }
    cout << "Test passed!" << endl << endl; 
#endif
    // ---

       cout << "Test " << ++numTests <<": test of SequenceReaderLocal" 
      << endl << endl;

          SequenceReaderLocal localReader(testReader,wordLength,cout);

          assert (    localReader.getNumSequencesInFile()
               ==  testReader.getNumSequencesInFile() );

            for ( int j(0) ; j < testReader.getNumSequencesInFile() ; j ++ ) 
        {
          wSingle.clear(); wMulti.clear(); 



        assert( testReader.getNextSequence(wSingle, wordLength) ==  
               localReader.getNextSequence(wMulti, wordLength) );
        testReader.getLastSequenceName(expected);
        testReader.getLastSequenceName(actual);
                cout << expected << " " << actual << endl;
        assert( expected == actual );
                 cout << localReader.getLastSequenceNumber() << " " 
                      << testReader.getLastSequenceNumber() << endl;
 
        assert(    localReader.getLastSequenceNumber()
                ==  testReader.getLastSequenceNumber() );
        assert ( wSingle == wMulti );
      } // ~for j

      assert( localReader.getNextSequence(wMulti, wordLength)== -1 );


        //    buffer << printWord(w,wordLength) << ends;
        //    actual=buffer.str();
        //  buffer.freeze(false); // release memory to write further stuff to buffer  
        //  buffer.seekp(0,ios::beg); // ensure next write is to start of buffer
 
    cout << "Test passed!" << endl << endl; 

    // ---

   cout << "Test " << ++numTests <<": test of SequenceReaderFilter" 
         << endl << endl;

   {

     SequenceReader* pFilteree = new SequenceReaderFastq("test_filter.fastq");
     SequenceReaderFilter filterer(pFilteree,"test_filter.fail");

     cout << "Check details match original SeqReader..." << endl;
     assert(filterer.getBitsPerSymbol()==pFilteree->getBitsPerSymbol());
     assert(filterer.getSourceDataType()==pFilteree->getSourceDataType());

     cout << "Check cycles through all sequences OK..." << endl;
     vector<WordSequence> seqs;
     vector<string> names;

     seqs.push_back(WordSequence());
     names.push_back(string());

     while (filterer.getNextSequence(seqs.back(),16)!=-1)
     {
       filterer.getLastSequenceName(names.back());
       cout << names.back() << endl << printWord(seqs.back(),16) << endl;
       seqs.push_back(WordSequence());
       names.push_back(string());
       
     }

     cout << "Read all sequences OK" << endl;

     cout << "Check no sequences have been \"lost\"..." << endl; 
     cout << pFilteree->getNumSequencesInFile()
          << " seqs in original file\n";
     cout << filterer.getNumSequencesInFile() 
          << " seqs in filtered file\n";
     cout << filterer.getNumFiltered() 
          << " filtered seqs\n";

     assert(    pFilteree->getNumSequencesInFile() 
             ==   filterer.getNumSequencesInFile() 
                + filterer.getNumFiltered() );

     
     cout 
       << "Check output of getSequence matches output of getNextSequence..."
       << endl;
     WordSequence seq;
     string name;

     for ( int i(1) ; i <= filterer.getNumSequencesInFile() ; i++ )
     {
       filterer.getSequence(seq, i, 16);
       filterer.getLastSequenceName(name);
       cout << name << endl << printWord(seq,16) << endl;
       assert(seq==seqs[i-1]);
       assert(name==names[i-1]);
     }

     cout << "Check rewind works..." << endl;
     filterer.rewind();

     assert (filterer.getNextSequence(seq,16)!=-1);

     assert (seq==seqs[0]);

     // Check won't accept dodgy sequence numbers
     assert( filterer.getSequence
             (seq, filterer.getNumSequencesInFile()+1, 16) 
             == -1 );
     assert( filterer.getSequence(seq, 0, 16) == -1 );



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

   // ---

   cout << "Test " << ++numTests 
        <<": test of SourceReader functionality of SequenceReaderFile" 
         << endl << endl;

   {
     SequenceReaderFasta testReader("test_extract.fasta");
     SequenceReaderMulti testReaderMulti(cout);
     testReaderMulti.addReader(testReader.clone());

     testReader.saveIndex("test_index1");
     testReaderMulti.saveIndex("test_index2");

     //     SourceReaderIndex testIndexReader1("test_index1");
     SourceReaderIndex testIndexReader("test_index2");


     string data=
       "tggtaaaaaaggttttcctgcatattaatatcatggtataggtgattgctgatggtatgtctttcttttactagatcatgcactgttcaggagcagaaaaatgtggctgttttttccaccacatattgtcctactactgtgtaccctgtctcagtaaatgccagtcagttacccaagctaaatgtccttgattcttccattctttcactacttccctggttcacattccagttggccaccaggtcatacttattttatcttaaaaattctctctctgaaatgcaaatcaaaactacgatgagataccatctcaaaccagttagtttgagatggcagttattaaaaagttcaaacataacagatgctgcgaggttgtggagaaaagaaaacacttatgcactgttggtgggagtgtaaattaattcaaccattgtggaaagtagtgtggtgatttctcaaagagctgaaagcagaactaccaatcaatccagccattccattactg-gtatatacccagaggaatataaattgttctatcataaagacacatgcataggtgtgtccattgtatcactactta-catagcaaagacatgaaatcaaccta-atgcccactggtgatagac";

     int seqSize=data.size(); 

     cout << seqSize << endl; 

     //     vector<char> v1, v2, v1m, v2m, v1i, v2i;
     char* v1; 
     char* v2; 
     char* v1m; 
     char* v2m; 
     char* v1i; 
     char* v2i; 

     string n1,n2;

     for (int j(1) ; j <= seqSize  ; j++)
     {
       cout << "j: " << j << endl;
       //      v1.clear(); v2.clear(); 
       //    v1m.clear(); v2m.clear();
       //   v1i.clear(); v2i.clear();

       testReader.extractSource(&v1, 1, j, seqSize );
       testReader.extractSource(&v2, 2, j, seqSize );

       testReaderMulti.extractSource(&v1m, 1, j, seqSize );
       testReaderMulti.extractSource(&v2m, 2, j, seqSize );

       testIndexReader.extractSource(&v1i, 1, j, seqSize );
       testIndexReader.extractSource(&v2i, 2, j, seqSize );

       //       assert(v1==v2);

       //       assert(v1==v1m);
       //   assert(v2==v2m);

       //       assert(v1==v1i);
       //     assert(v2==v2i);

       for ( int (i(j)) ; i <=seqSize ; i++ ) assert(v1[i-j]==data[i-1]);
       for ( int (i(j)) ; i <=seqSize ; i++ ) assert(v2[i-j]==data[i-1]);
       for ( int (i(j)) ; i <=seqSize ; i++ ) assert(v1m[i-j]==data[i-1]);
       for ( int (i(j)) ; i <=seqSize ; i++ ) assert(v2m[i-j]==data[i-1]);
       for ( int (i(j)) ; i <=seqSize ; i++ ) assert(v1i[i-j]==data[i-1]);
       for ( int (i(j)) ; i <=seqSize ; i++ ) assert(v2i[i-j]==data[i-1]);

     }
     cout << "got through  first bit" << endl;

     for (int j(1) ; j <= seqSize  ; j++)
     {
       //       v1.clear(); v2.clear(); 
       //     v1m.clear(); v2m.clear();
       //     v1i.clear(); v2i.clear();

       testReader.extractSource(&v1, 1, 1, j );
       testReader.extractSource(&v2, 2, 1, j );

       testReaderMulti.extractSource(&v1m, 1, 1, j );
       testReaderMulti.extractSource(&v2m, 2, 1, j );

       testIndexReader.extractSource(&v1i, 1, 1, j );
       testIndexReader.extractSource(&v2i, 2, 1, j );

       //       assert(v1==v2);

       //       assert(v1==v1m);
       //   assert(v2==v2m);

       //       assert(v1==v1i);
       //   assert(v2==v2i);

       for ( int i(1) ; i < j ; i++ ) assert(v1[i]==data[i]);
       for ( int i(1) ; i < j ; i++ ) assert(v2[i]==data[i]);
       for ( int i(1) ; i < j ; i++ ) assert(v1m[i]==data[i]);
       for ( int i(1) ; i < j ; i++ ) assert(v2m[i]==data[i]);
       for ( int i(1) ; i < j ; i++ ) assert(v1i[i]==data[i]);
       for ( int i(1) ; i < j ; i++ ) assert(v2i[i]==data[i]);

     }
     cout << "got through second bit" << endl;



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

   // ---

   cout << "Test " << ++numTests 
        <<": test that SourceReader functionality works with getNextSequence" 
         << endl << endl;

   {
     string data=
       "tggtaaaaaaggttttcctgcatattaatatcatggtataggtgattgctgatggtatgtctttcttttactagatcatgcactgttcaggagcagaaaaatgtggctgttttttccaccacatattgtcctactactgtgtaccctgtctcagtaaatgccagtcagttacccaagctaaatgtccttgattcttccattctttcactacttccctggttcacattccagttggccaccaggtcatacttattttatcttaaaaattctctctctgaaatgcaaatcaaaactacgatgagataccatctcaaaccagttagtttgagatggcagttattaaaaagttcaaacataacagatgctgcgaggttgtggagaaaagaaaacacttatgcactgttggtgggagtgtaaattaattcaaccattgtggaaagtagtgtggtgatttctcaaagagctgaaagcagaactaccaatcaatccagccattccattactg-gtatatacccagaggaatataaattgttctatcataaagacacatgcataggtgtgtccattgtatcactactta-catagcaaagacatgaaatcaaccta-atgcccactggtgatagac";

     SequenceReaderFasta testReader("test_extract.fasta");
     SequenceReaderMulti testReaderMulti(cout);
     testReaderMulti.addReader(testReader.clone());
     SourceReaderIndex testIndexReader("test_index2");

     WordSequence w1, w2;
     //     vector<char> v1, v2;
     char* v1; char* v2;

     int numSeqs(4);
     int numChars(100);

     string n1,n2;

     for (int i(0);i<numSeqs; i++)
     {
       cout << "i: " << i << endl;
       w1.clear(); w2.clear();

       assert(testReader.getLastSequenceNumber()==i);
       assert(testReaderMulti.getLastSequenceNumber()==i);

       assert(testReader.findSequence(i+1)==true);
       assert(testReaderMulti.findSequence(i+1)==true);

       assert(testReader.getLastSequenceNumber()==i);
       assert(testReaderMulti.getLastSequenceNumber()==i);

     
       assert(testReader.getNextSequence( w1, 10 )!=-1);
       assert(testReaderMulti.getNextSequence( w2, 10 )!=-1);
       assert(w1==w2);
    
       assert(testReader.getLastSequenceNumber()==i+1);
       assert(testReaderMulti.getLastSequenceNumber()==i+1);
       //       cout << "fred: " << testReader.getLastSequenceNumber()
       // << " " << testReaderMulti.getLastSequenceNumber() << endl;
       testReader.getLastSequenceName(n1);
       testReaderMulti.getLastSequenceName(n2);
       assert(n1==n2);
       cout << n1 << " " << n1.size() << endl;
       cout << testIndexReader.extractName(i+1) << endl;
       n2=(string)testIndexReader.extractName(i+1);
       cout << n2 << " " << n2.size() << endl;
       assert(n1==n2);

       for ( int j(1); j <= numSeqs ; j++ )
       {
         cout << "j: " << j << endl;
         //      v1.clear(); v2.clear();
         testReader.extractSource(&v1, j, 1, numChars );
         testReaderMulti.extractSource(&v2, j, 1, numChars );
         //      assert(v1==v2);

         //       cout << "fred2: " << testReader.getLastSequenceNumber()
         // << " " << testReaderMulti.getLastSequenceNumber() << endl;

         assert(testReader.getLastSequenceNumber()==i+1);
         assert(testReaderMulti.getLastSequenceNumber()==i+1);

         for ( int k(0) ; k < numChars ; k++ ) assert(v1[k]==data[k]);
         for ( int k(0) ; k < numChars ; k++ ) assert(v2[k]==data[k]);


       } // ~for j

     } // ~for i







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

   // ---





    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 testSequenceReaderFasta




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

// End of file testSequenceReaderFasta.cpp

---