Global/SSAHAMain.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 2004

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

// Description:

// Includes:

#include "SSAHAMain.h"
#include "SequenceEncoder.h"
#include "SequenceReaderMulti.h"
#include "SequenceReaderFasta.h"
#include "SequenceReaderFastq.h"
//#include "SequenceReaderCodon.h"
#include "SequenceReaderFilter.h"
#include "HashTable.h"
#include "HashTablePacked.h"
#include "HashTableTranslated.h"
#include "MatchStore.h"
#include "MatchStoreUngapped.h"
#include "MatchStoreGapped.h"
#include "MatchAligner.h"
#include "TimeStamp.h"
#include <dirent.h>
#include <fstream>
#include <string>
#include <assert.h>
#include <exception>
#include <cmath>

MachineInfo* MachineInfo::info_ = 0;
NullBuffer nullBuffer;
ofstream*  pLogStream = NULL;

int defaultWordLengthDNA = 12;
int defaultWordLengthProtein = 4;



// ### Function Definitions ###

//typedef MatchTaskSort
//<CombineSort<SortByMatchSize,SortBySubjectName,SortByQueryStart> >
//SortByMatchType;

typedef MatchTaskSort
<SortByMatch>
SortByMatchType;





int main( int numArgs, char* args[] )
{

  try
  {

    cerr << 
    "Welcome to SSAHA: Sequence Search and Alignment by Hashing Algorithm\n\
Copyright (C) 2004 by Genome Research Limited\n\
This software is released under the terms of version 2 of the GNU General\n\
Public Licence, as published by the Free Software Foundation.\n\
This is SSAHA Version 3.2, released 1st March 2004.\n\n";

  Timer timeStamp;

  string thisArg;
  if ( numArgs < 3 )
  {
    // -h, -help or no argument prints help message
    // anything else throws an exception
    if ( numArgs == 2 )
    {
      thisArg = (string) args[1];
      if ( ( thisArg != "-h" ) && ( thisArg != "-help" ) )
      {
        cerr << "Error: invalid argument"
             << " (" << args[1] << ")" 
             << " - use \'" << args[0] << " -help\' for help.\n";
        throw SSAHAException("Invalid command line argument.");
      } // ~if
    } // ~if


    int i(0); 
    do 
    {
      cout << helpText[i] << "\n";
    } while (helpText[++i][0]!='@');

    return(0);



  } // ~if

  QueryParameterStruct queryParams( defaultParams );

  CommandLineArg::parseCommandLine( numArgs, args, queryParams );

  //  cerr << timeStamp;

  processQuery( queryParams );

  //  cerr << timeStamp << "Job finished.\n";

  if (pLogStream != NULL) delete pLogStream;

  return (0);

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

} // ~main

// Name:      parseCommandLine
// Arguments: command line arguments (in), QueryParameterStruct& (out)
// Go through the command line arguments: if an argument is a recognized
// keyword, take the following argument and place it in the appropriate
// part of queryParams. Does not make any attempt to check if paramter values
// are sensible.
void CommandLineArg::parseCommandLine
( int numArgs, char* args[], QueryParameterStruct& queryParams )
{

  int firstArg;

  // Add mandatory parameters
  if ( (args[2])[0] == '-' )
  {
    // then only one file name given: assume this is a subject file
    // to be hashed
    firstArg = 2;
    queryParams.subjectName = (string) args[1];
    queryParams.runQuery = false; 
  }
  else
  {
    // two file names given: assume first is query, second is subject
    firstArg = 3;
    queryParams.queryName   = (string) args[1];
    queryParams.subjectName = (string) args[2];
  }

  // Parse options
  vector<CommandLineArg*> validArgs;
  validArgs.push_back
  ( new CommandLineArgInt( "-wordLength", "-wl", queryParams.wordLength ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-stepLength", "-sl", queryParams.stepLength ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-maxStore", "-ms", queryParams.maxStore ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-minPrint", "-mp", queryParams.minPrint ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-maxGap", "-mg", queryParams.maxGap ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-maxInsert", "-mi", queryParams.maxInsert ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-numRepeats", "-nr", queryParams.numRepeats ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-substituteWords", "-sw", 
                           queryParams.substituteThreshold ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-bandExtension", "-be", 
                           queryParams.bandExtension ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-sortMatches", "-sm", queryParams.sortMatches ) );
  validArgs.push_back
  ( new CommandLineArgInt
    ( "-doAlignment", "-da", queryParams.doAlignment ));
  validArgs.push_back
  ( new CommandLineArgInt( "-queryStart", "-qs", queryParams.queryStart ) );
  validArgs.push_back
  ( new CommandLineArgInt( "-queryEnd", "-qe", queryParams.queryEnd ) );
  validArgs.push_back
  ( new CommandLineArgBool( "-hashStats", "-hs", queryParams.printHashStats ));
  validArgs.push_back
  ( new CommandLineArgBool( "-packHits", "-ph", queryParams.packHits ));
  validArgs.push_back
  ( new CommandLineArgBool
    ( "-reverseQuery", "-rq", queryParams.reverseQuery ));  
  validArgs.push_back
  ( new CommandLineArgBool
    ( "-parserFriendly", "-pf", queryParams.parserFriendly ) );
  validArgs.push_back
  ( new CommandLineArgString( "-queryType", "-qt", queryParams.queryType ) );
  validArgs.push_back
  ( new CommandLineArgString( "-subjectType", "-st", queryParams.subjectType));
  validArgs.push_back
  ( new CommandLineArgString( "-queryFormat", "-qf", queryParams.queryFormat));
  validArgs.push_back
  ( new CommandLineArgString
    ( "-subjectFormat", "-sf", queryParams.subjectFormat) );
  //  validArgs.push_back
  //  ( new CommandLineArgString( "-reportMode", "-rm", queryParams.reportMode ) );
  validArgs.push_back
  ( new CommandLineArgString( "-queryReplace", "-qr", queryParams.queryReplace ) );
  validArgs.push_back
  ( new CommandLineArgString( "-subjectReplace", "-sr", queryParams.subjectReplace ) );
  validArgs.push_back
  ( new CommandLineArgString( "-logMode", "-lm", queryParams.logMode ) );
  validArgs.push_back
  ( new CommandLineArgString( "-saveName", "-sn", queryParams.saveName ) );

  string thisArg;
  vector<CommandLineArg*>::iterator pArg = validArgs.end();
  for ( int i(firstArg); i < numArgs ; i++ )
  {
    thisArg = (string) args[i];
    if (pArg != validArgs.end()) 
    // ... then an instance of CommandLineArg is waiting to parse the current
    // command line arg
    { 
      (*pArg)->addValue( thisArg ); 
      pArg = validArgs.end(); 
    } // ~if
    else
    // ... trawl through the list of command line args and see if any of 'em
    // recognize the current command line arg
    { 
      for (pArg = validArgs.begin() ; pArg != validArgs.end(); ++pArg )
      {
        if ( (*pArg)->isThisMe(thisArg) ) break; // *pArg recognizes the arg
      } // ~for      
      if ( pArg == validArgs.end() )
      // ... then the argument has not been recognized - complain!
      {
        cerr << "Error: command line argument " << thisArg << " not found.\n";
        exit(1);
      } // ~if

      // if pArg is an instance of type CommandLineArgBool, then there is
      // no associated value, so don't try to read one
      if ( dynamic_cast<CommandLineArgBool*>(*pArg)) pArg=validArgs.end();
    } // ~else    
  } // ~for

  // Deallocate memory
  for ( vector<CommandLineArg*>::iterator i = validArgs.begin(); 
        i != validArgs.end(); ++i ) delete *i;


} // ~parseCommandLine

// Name:      processQuery
// Arguments: QueryParameterStruct& (in)
// 1. Takes queryParams produced by parseCommandLine
// 2. Check the values therein are sensible.
// 3. Set up a SequenceReader for the query data
// 4. Either load in a hash table or create one from the subject data
// 5. Pass the hash table and query sequence reader to the QueryManager
// 6. Run the query!
void processQuery( QueryParameterStruct& queryParams )
{

  // Send log information to appropriate destination
  if ( queryParams.logMode == "cout" )
  {
    cerr.rdbuf(cout.rdbuf());
  } // ~if
  else if ( queryParams.logMode == "null" )
  {
    cerr.rdbuf(&nullBuffer);
  } // ~else if
  else if ( queryParams.logMode != "cerr" )
  {
    pLogStream = new ofstream(queryParams.logMode.c_str(), ios::ate|ios::app);
    if (pLogStream->fail())
    {
      cerr << "could not open log stream file " << queryParams.logMode
           << ".\n";
      throw SSAHAException("could not open log file.\n");
    } // ~if
    cerr.rdbuf(pLogStream->rdbuf());
  } // ~else if 

  cerr << "Run started: " << getTimeNow()
       << "Input query parameters:\n" << queryParams;

  // Obtain and print machine info
  MachineInfo* info = MachineInfo::getInfo();
  cerr << *info;

  // Check queryParams values are sensible

  if ( queryParams.minPrint <= 0 )
  {
    cerr << "Error: value for minPrint (" << queryParams.minPrint 
         << ") must be greater than zero.\n";
    throw SSAHAException("Invalid value for minPrint");
  } // ~if

  if ( queryParams.maxStore < 0 )
  {
    cerr << "Error: value for maxStore (" << queryParams.maxStore 
         << ") must be greater than zero.\n";
    throw SSAHAException("Invalid value for maxStore");
  } // ~if
  if ( queryParams.substituteThreshold < 0 )
  {
    cerr << "Error: value for substituteThreshold (" 
         << queryParams.substituteThreshold
         << ") must be greater than zero.\n";
    throw SSAHAException("Invalid value for substituteThreshold");
  } // ~if


  if ( queryParams.doAlignment != 0 )
  {
    if( queryParams.doAlignment < 20 )
    {
      cerr << "Error: value for doAlignment (" << queryParams.doAlignment 
           << ") must be at least 20.\n";
      throw SSAHAException("Invalid value for doAlignment");
    } // ~if
#ifdef XXX
    if (    ( queryParams.reportMode != "replaceA" )
         && ( queryParams.reportMode != "tag" )   )
    {
      cerr << "Info: graphical alignments can't be produced with reportMode\n"
           << "set to \"" << queryParams.reportMode 
           << "\" , setting reportmode to \"tag\" and proceeding.\n";
      queryParams.reportMode="tag";
    } // ~if   
#endif
    if (    ( queryParams.queryReplace == "ignore" )
         || ( queryParams.queryReplace == "report" )   )
    {
      cerr << "Info: graphical alignments can't be produced with queryReplace\n"
           << "set to \"" << queryParams.queryReplace 
           << "\" , resetting to default and proceeding.\n";
      queryParams.queryReplace="default";
    } // ~if   

    if (    ( queryParams.subjectReplace == "ignore" )
         || ( queryParams.subjectReplace == "report" )   )
    {
      cerr << "Info: graphical alignments can't be produced with subjectReplace\n"
           << "set to \"" << queryParams.subjectReplace
           << "\" , setting to \"tag\" and proceeding.\n";
      queryParams.subjectReplace="tag";
    } // ~if   


  } // ~if

  // If no step length specified, set it to the word length
  if (    ( queryParams.stepLength <= 0 )
       || ( queryParams.stepLength > queryParams.wordLength )   )
  { 
    cerr << "Info: setting stepLength to " << queryParams.wordLength << ".\n";
    queryParams.stepLength = queryParams.wordLength;
  } // ~if

  // Check that queryStart and queryEnd are sensible
  if (queryParams.queryStart <= 0 )
  {
      cerr << "Error: first requested sequence  (" 
      << queryParams.queryStart
      << ") must be greater than or equal to 1.\n";
      throw SSAHAException("Invalid query start value!");
  } // ~if

  if (    ( queryParams.queryEnd != - 1 ) 
       && ( queryParams.queryEnd < queryParams.queryStart ) )
  {
      cerr << "Error: last requested sequence  (" 
      << queryParams.queryEnd
      << ") should be greater than or equal to first ("
      << queryParams.queryStart << ").\n";
      throw SSAHAException("Invalid query start and end values!");
  } // ~if

  // Set up mode for sequence readers

#ifdef XXXX
  SequenceReaderMode* pMode;

  if ( queryParams.reportMode == "tag" )
  {
    pMode = new SequenceReaderModeFlagReplace('A');
  }
  else if ( queryParams.reportMode == "ignore" )
  {
    pMode = new SequenceReaderModeIgnore();
  }
  else if ( queryParams.reportMode == "report" )
  {
    pMode = new SequenceReaderModeReport();
  }
  else if ( queryParams.reportMode == "replaceA" )
  {
    pMode = new SequenceReaderModeReplace('A');
  }
  else if ( queryParams.reportMode == "replaceG" )
  {
    pMode = new SequenceReaderModeReplace('G');
  }
  else if ( queryParams.reportMode == "replaceC" )
  {
    pMode = new SequenceReaderModeReplace('C');
  }
  else if ( queryParams.reportMode == "replaceT" )
  {
    pMode = new SequenceReaderModeReplace('T');
  }
  else if ( queryParams.reportMode == "rrepA" )
  {
    pMode = new SequenceReaderModeReportReplace('A');
  }
  else if ( queryParams.reportMode == "rrepG" )
  {
    pMode = new SequenceReaderModeReportReplace('G');
  }
  else if ( queryParams.reportMode == "rrepC" )
  {
    pMode = new SequenceReaderModeReportReplace('C');
  }
  else if ( queryParams.reportMode == "rrepT" )
  {
    pMode = new SequenceReaderModeReportReplace('T');
  }
  else 
  {
    cerr << "Error: value for reportMode (" << queryParams.reportMode 
         << ") not recognised.\n";
    throw SSAHAException("Invalid value for reportMode");
  } // ~else
#endif

  // trap foolish command line argument combinations
  if (    ( queryParams.subjectFormat == "hash" )
       && ( queryParams.saveName != "" ) )
  {
    cerr << "Error: trying to create hash table from hash table!\n";
    throw SSAHAException("Invalid value for subject format");
  } // ~if

  if ( queryParams.runQuery == false )
  {
    cerr << 
    "Info: no query name given, assuming hash table generation only \
is required.\n";
    if (    ( queryParams.saveName == "" )
         && ( queryParams.printHashStats==false )   )
    {
      cerr << "Error: no saveName for hash table has been given.\n";
      throw SSAHAException("No saveName for hash table");
    } // ~if
  } // ~if

  // Set up subject data ...
  //  Allocator<PositionInHitList>* pArrayAllocator;
  //  cerr << "Info: pointer array will be allocated ";

  //  if (queryParams.mapArray==true)
  //  {
  //   cerr << "by memory mapping.\n";
  //    pArrayAllocator = new AllocatorMapped<PositionInHitList>;
  // }
  //  else
  //  {
  //   cerr << "from local memory.\n";
  //   pArrayAllocator = new AllocatorLocal<PositionInHitList>;
  //  }

  //  Allocator<PositionInDatabase>* pHitListAllocator;
  //  Allocator<PositionPacked>* pHitListAllocator;
  //  cerr << "Info: hit list will be allocated ";

  //  if (queryParams.mapHits==true)
  //  {
  //   cerr << "by memory mapping.\n";
  //   pHitListAllocator = new AllocatorMapped<PositionInDatabase>;
  //  }
  //  else
  //  {
  //   cerr << "from local memory.\n";
    //    pHitListAllocator = new AllocatorLocal<PositionInDatabase>;
  //   pHitListAllocator = new AllocatorLocal<PositionPacked>;
  //  }

  //  HashTable hashTable

  HashTableFactory creator(cerr);

  HashTableGeneric* pHashTable(NULL);
  SourceReader* pSubjectSource(NULL);

  if ( queryParams.subjectFormat == "hash" ) 
  { 
    SourceReaderIndex* pIndexReader;
    // loading in a previously created hash table

    {
      ifstream filesFile( (queryParams.subjectName+(string)".files").c_str() );
      ifstream indexFile( (queryParams.subjectName+(string)".index").c_str() );
      if ((!filesFile.fail())&&(!indexFile.fail()))
      {
        pIndexReader = new SourceReaderIndex(queryParams.subjectName);
        cerr << "Info: found a sequence source index for this hash table.\n";
      } // ~if  
      else
      {
        pIndexReader = NULL;
        cerr << 
          "Info: did not find a sequence source index for this hash table.\n";
            
      } // ~else
     filesFile.close();
     indexFile.close();
   } // ~scope of ifstreams

    pHashTable = creator.loadHashTable
      ( queryParams.subjectName, pIndexReader );
    pHashTable->setMaxNumHits( queryParams.maxStore );

    pSubjectSource=pIndexReader; // this ensures it gets deleted at the end

    // Determine and display type of loaded in hash table
    // subjectType still needs to be properly set up as it is
    // used to determind which MatchTaskAlign is needed
    if ( pHashTable->getBitsPerSymbol() == gResidueBits )
    {
      if ( pHashTable->getSourceDataType() == gProteinData )
      {
        cerr << "Hash table contains protein data." << endl;
        queryParams.subjectType = "protein";
      } // ~if
      else if ( pHashTable->getSourceDataType() == gDNAData )
      {
        cerr << "Hash table contains translated DNA data." << endl;
        queryParams.subjectType = "codon";
      } // ~else if
      else 
      {
        cerr << "Warning: could not determine type of hash table, "
             << "attempting to continue..." << endl;
      } // ~else
    } // ~if
    else if (     ( pHashTable->getBitsPerSymbol() == gBaseBits )
               && ( pHashTable->getSourceDataType() == gDNAData ) ) 
    {
        cerr << "Hash table contains untranslated DNA data." << endl;
        queryParams.subjectType = "DNA";
    } // ~if
    else
    {
        cerr << "Warning: could not determine type of hash table, "
             << "attempting to continue..." << endl;
    } // ~else


    // any word length / step length input is irrelevant as
    // we have to use the values the table was created with.
    if ( pHashTable->getWordLength() != queryParams.wordLength )
    {
      cerr << "Info: hash table was created using hash word length of "
           << pHashTable->getWordLength() 
           <<" bases,\nproceeding using this value.\n";
      queryParams.wordLength = pHashTable->getWordLength();
    } // ~if

    if ( pHashTable->getStepLength() != queryParams.stepLength )
    {
      cerr << "Info: hash table was created using interval of "
           << pHashTable->getStepLength() 
           <<" between bases,\nproceeding using this value.\n";
      queryParams.stepLength = pHashTable->getStepLength();
    } // ~if


  } // ~if
  else 
  {
    // have to create our own hash table

    if (    ( queryParams.subjectType == "DNA" ) 
         && ( queryParams.queryType == "protein" ) )
    {
      cerr << "Warning: can't run protein query against DNA database, setting "
           << "subject database type to \"codon\" instead." << endl;
      queryParams.subjectType = "codon";
    } // ~if

    if ( queryParams.subjectType == "DNA" )
    {
      pHashTable = new HashTablePacked( cerr, queryParams.saveName );
      if (    ( queryParams.wordLength <= 0 )
              || ( queryParams.wordLength*gBaseBits > ( 8*sizeof(Word)) -1 ) )
      {
        cerr << "Warning: word length (" << queryParams.wordLength 
             << ") outside valid range (0 to " 
             << (8*sizeof(Word)-1)/gBaseBits << "),\n"
             << "using default word length of " 
             << defaultWordLengthDNA << " instead.\n";
        queryParams.wordLength = defaultWordLengthDNA;
        //      throw SSAHAException("Invalid value for wordLength.");
      
      } // ~if
      
    } // ~if
    else 
    {
      if ( queryParams.subjectType == "protein" )
      {
        
        pHashTable = new HashTablePackedProtein( cerr, queryParams.saveName );
        //      pHashTable = new HashTableFred( cerr, queryParams.saveName );
      } // ~if
      else if ( queryParams.subjectType == "codon" )
      {
        pHashTable = new HashTableTranslated( cerr, queryParams.saveName );
      } // ~else if
      else 
      {
        cerr << "Error: value for subjectType (" << queryParams.subjectType
             << ") not recognized.\n";
        throw SSAHAException("Invalid value for subjectType.\n");
      } // ~if

      if (    ( queryParams.wordLength <= 0 )
              || ( queryParams.wordLength*gResidueBits > ( 8*sizeof(Word)) -1 ) )
      {
        cerr << "Warning: word length (" << queryParams.wordLength 
             << ") outside valid range (0 to " 
             << (8*sizeof(Word)-1)/gResidueBits << "),\n"
             << "using default word length of " 
             << defaultWordLengthProtein << " instead.\n";
        queryParams.wordLength = defaultWordLengthProtein;
        //      throw SSAHAException("Invalid value for wordLength.");
      } // ~if
     } // ~else

    // If no step length specified, set it to the word length
    if (    ( queryParams.stepLength <= 0 )
         || ( queryParams.stepLength > queryParams.wordLength )   )
    { 
      cerr << "Info: setting stepLength to " << queryParams.wordLength << ".\n";
      queryParams.stepLength = queryParams.wordLength;
    } // ~if
    

    SequenceReaderMode* pSubjectMode(NULL);

    if ( queryParams.subjectReplace == "tag" )
    {
      pSubjectMode = new SequenceReaderModeFlagReplace('A');
    } // ~if
    else if ( queryParams.subjectReplace == "ignore" )
    {
      pSubjectMode = new SequenceReaderModeIgnore();
    } // ~else if
    else if ( queryParams.subjectReplace == "report" )
    {
      pSubjectMode = new SequenceReaderModeReport();
    } // ~else if
    else if ( queryParams.subjectReplace.size()==1)
    {
      unsigned char replaceChar(*queryParams.subjectReplace.c_str());
      if ((queryParams.subjectType=="protein")&&(ttProtein[replaceChar]==nv))
      {
        cerr << "Error: value for subjectReplace (" << queryParams.subjectReplace 
             << ") must be a valid IUPAC amino acid code.\n";
        throw SSAHAException("Invalid value for subjectReplace");
      } // ~if
      else if (ttDNA[replaceChar]==nv)
      {
        cerr << "Error: value for subjectReplace (" << queryParams.subjectReplace 
             << ") must be A,G,C or T.\n";
        throw SSAHAException("Invalid value for subjectReplace");
      } // ~else if
      pSubjectMode = new SequenceReaderModeReplace(replaceChar);
    } // ~else if
  else 
  {
    cerr << "Error: value for subjectReplace (" << queryParams.subjectReplace 
         << ") not recognised.\n";
    throw SSAHAException("Invalid value for subjectReplace");
  } // ~else


    SequenceEncoder* pSubjectEncoder;

    if ( queryParams.subjectType == "protein" )
      pSubjectEncoder = new SequenceEncoderProtein(5,cerr);
    else
      pSubjectEncoder = new SequenceEncoderDNA(12,cerr); 

    SequenceReader* pSubject =
    findSequenceReader
    ( queryParams.subjectName, queryParams.subjectFormat, pSubjectEncoder ); 
    assert (pSubject !=NULL);

    pSubject->changeMode( pSubjectMode );
    delete pSubjectMode; // cloned, so no longer needed

    creator.createHashTable
    (  *pHashTable,
       *pSubject, 
       queryParams.wordLength, 
       queryParams.maxStore, 
       queryParams.stepLength  );

    // Save the hash table to a file, if required
    if (    ( queryParams.saveName != "" )
         && ( queryParams.subjectFormat != "hash" )  )
    { 
      //      hashTable.saveHashTable(queryParams.saveName);
      creator.saveHashTable( *pHashTable );
      pSubject->saveIndex( queryParams.saveName );
    } // ~if

    if (queryParams.doAlignment != 0 )
    {
      pSubjectSource = (SourceReader*) pSubject;
    } // ~if
    else
    {
      delete pSubjectEncoder;
      delete pSubject;
    } // ~else
  } // ~else

  assert (pHashTable!=NULL);
  assert (pHashTable->isInitialized());

  double expectedNumHits = pHashTable->getTotalNumWords();

  if (    (pHashTable->getHitListFormat()==gTranslated)
       || (pHashTable->getHitListFormat()==g32BitPackedProtein) )
  {
    expectedNumHits /= 
      pow((double)gNumCodonEncodings,(int)pHashTable->getWordLength());
  } // ~if
  else
  {
    assert(pHashTable->getBitsPerSymbol()==gBaseBits);
    expectedNumHits /= 1<<(2*pHashTable->getWordLength());
  } // ~else

  //  expectedNumHits /= (pHashTable->getBitsPerSymbol()==gBaseBits)
  //   ? 1<<(2*pHashTable->getWordLength())
  //   : pow((double)gNumCodonEncodings,(int)pHashTable->getWordLength());

  cerr << "Info: would expect " << expectedNumHits 
         << " hits per word for a random database of this size." << endl;

  queryParams.maxStore=1+(int)(expectedNumHits*queryParams.maxStore);

  cerr << "Info: will ignore hits on words that occur more than " 
       << queryParams.maxStore << " times in the database." << endl;

  pHashTable->setMaxNumHits
  (
    queryParams.maxStore
  );


  if ( queryParams.printHashStats == true )
  {
    pHashTable->printHashStats();
  } // ~if

  if ( queryParams.runQuery == false )
  {
    cerr << "Info: no query name given, assuming no query required.\n"; 
    delete pHashTable;
    return;
  } // ~if

  // Check that numRepeats is sensible
  // Checked here and not earlier because word length may change from
  // the value specified in the command line arguments if the hash
  // table has been loaded in from disk.
  if ( queryParams.numRepeats < 0 )
  {
    cerr << "Error: max size of repeated motif to be masked (" 
         << queryParams.numRepeats
         << ") outside valid range (0 to " 
         << queryParams.wordLength << ").\n";
    throw SSAHAException("Invalid value for numRepeats!");
  } // ~if
  else if ( queryParams.numRepeats > queryParams.wordLength ) 
  {
    cerr << "Warning: max size of repeated motif to be masked (" 
         << queryParams.numRepeats
         << ") outside valid range (0 to " 
         << queryParams.wordLength << "),\nusing "
         << queryParams.wordLength << " instead.\n";
    queryParams.numRepeats = queryParams.wordLength;
  } // ~if

  if (    ( queryParams.queryType != "DNA" )
          && ( queryParams.queryType != "protein" )
          && ( queryParams.queryType != "codon" ) )
  {
    cerr << "Error: value for queryType (" << queryParams.queryType
         << ") not recognized.\n";
    throw SSAHAException("Invalid value for queryType.\n");
  } // ~if

  if (queryParams.substituteThreshold > 0 )
  {
    pHashTable->setSubstituteThreshold( queryParams.substituteThreshold );
  } // ~if

  SequenceEncoder* pQueryEncoder;
  
if ( queryParams.queryType == "protein" )
{
  if (    (dynamic_cast<HashTablePackedProtein*>(pHashTable)==NULL)
          && (dynamic_cast<HashTableTranslated*>(pHashTable)==NULL) )
  {
    cerr << "Error: protein query can't run against DNA database!\n";
    throw SSAHAException("Data type mismatch between query and subject");
  } // ~if
  pQueryEncoder = new SequenceEncoderProtein(5);
  if (queryParams.queryReplace == "default") queryParams.queryReplace="X";
} // ~if
else
{
  pQueryEncoder = new SequenceEncoderDNA(12); 
  if (queryParams.queryReplace == "default") queryParams.queryReplace="A";
} // ~else

assert (pQueryEncoder != NULL);


SequenceReaderMode* pQueryMode(NULL);

if ( queryParams.queryReplace == "tag" )
{
  pQueryMode = new SequenceReaderModeFlagReplace('A');
} // ~if
else if ( queryParams.queryReplace == "ignore" )
{
  pQueryMode = new SequenceReaderModeIgnore();
} // ~else if
else if ( queryParams.queryReplace == "report" )
{
  pQueryMode = new SequenceReaderModeReport();
} // ~else if
else if ( queryParams.queryReplace.size()==1)
{
  unsigned char replaceChar(*queryParams.queryReplace.c_str());
  if (queryParams.queryType=="protein")
  {
    if (ttProtein[replaceChar]==nv)
    {
      cerr << "Error: value for queryReplace (" << queryParams.queryReplace 
           << ") must be a valid IUPAC amino acid code.\n";
      throw SSAHAException("Invalid value for queryReplace");
    } // ~if
  } // ~if
  else if (ttDNA[replaceChar]==nv)
  {
    cerr << "Error: value for queryReplace (" << queryParams.queryReplace 
         << ") must be A,G,C or T.\n";
    throw SSAHAException("Invalid value for queryReplace");
  } // ~else if
  pQueryMode = new SequenceReaderModeReplace(replaceChar);
} // ~else if
else 
{
  cerr << "Error: value for queryReplace (" << queryParams.queryReplace 
       << ") not recognised.\n";
    throw SSAHAException("Invalid value for queryReplace");
} // ~else



  // Set up SequenceReader for query data
  SequenceReader* pQuery = 
  findSequenceReader
    ( queryParams.queryName, queryParams.queryFormat, pQueryEncoder );
  assert(pQuery!=NULL);  

  pQuery->changeMode( pQueryMode );
delete pQueryMode;

  // Set up query manager
  QueryManager queryManager( *pQuery,*pHashTable,cerr);

  // Set up MatchTask

  MatchTask* pTask(NULL);
  MatchTask* pPrintTask(NULL);
  //  MatchTaskPrint printNormal;
  //  MatchTaskPrintTabbed printTabbed;
  SortByMatchType sorter(queryParams.sortMatches,0.25);
  
  if (queryParams.doAlignment != 0 )
  {
    MatchAligner* pAligner;

    if (queryParams.queryType == "DNA")
    {
      if (queryParams.subjectType == "DNA")
      {
        pAligner = new MatchAlignerDNA
          ( queryParams.doAlignment, queryParams.bandExtension ); 
      } // ~if
      else if (queryParams.subjectType == "protein")
      {
        pAligner = new MatchAlignerTranslatedProtein
          ( false, queryParams.doAlignment, queryParams.bandExtension );   
      } // ~else if
      else // (queryParams.subjectType == "codon")
      {
        pAligner = new MatchAlignerTranslatedDNA
          ( queryParams.doAlignment, queryParams.bandExtension );   
      } // ~else

    } // ~if (queryParams.queryType == "DNA")
    else // query is protein
    {  
      if (queryParams.subjectType == "protein")
      {
        pAligner = new MatchAlignerProtein
          ( queryParams.doAlignment, queryParams.bandExtension ); 
      } // ~if
      else // subjectType is "DNA" or "codon"
      {
        pAligner = new MatchAlignerTranslatedProtein
          ( true, queryParams.doAlignment, queryParams.bandExtension );   
      } // ~else

    } // ~else

    assert( pAligner!=false);

    // ownership of *pAligner passes to *pPrintTask
    pPrintTask = new MatchTaskAlign( *pQuery, *pSubjectSource, pAligner,
                                     true, true );

  } // ~if
  else
  {

    if (queryParams.parserFriendly==true)
      {
        if (queryParams.reverseQuery==true) 
          pPrintTask = new MatchTaskPrintTabbedReverse(cout);
        else
          pPrintTask = new MatchTaskPrintTabbed(cout);
      } // ~if
    else
      {
        if (queryParams.reverseQuery==true) 
          pPrintTask = new MatchTaskPrintReverse(cout);
        else
          pPrintTask = new MatchTaskPrint(cout);
      } // ~else

  } // ~else

  assert(pPrintTask!=NULL);

  if (queryParams.sortMatches==0)
    pTask=pPrintTask;
  else
    pTask=new CombineTaskVirtual(sorter,*pPrintTask);

  assert (pTask!=NULL);

  if ((queryParams.maxGap == -1)&&(queryParams.maxInsert==0))
  {

    // This no longer works with the new codon stuff TC 11.12.1
    //    MatchAlgorithmUngapped algorithm
    //    ( queryParams.minPrint, queryParams.numRepeats );

    MatchAlgorithmGapped algorithm
    ( 0, 0, queryParams.minPrint, 
      queryParams.numRepeats );    

    queryManager.doQuery
    ( algorithm, *pTask, queryParams.queryStart, queryParams.queryEnd );


  } // ~if
  else
  {

    if (queryParams.maxGap==-1) queryParams.maxGap=0;

    MatchAlgorithmGapped algorithm
    ( queryParams.maxGap, queryParams.maxInsert, queryParams.minPrint, 
      queryParams.numRepeats );    

    queryManager.doQuery
    ( algorithm, *pTask, queryParams.queryStart, queryParams.queryEnd );
    

  } // ~else

  // Don't need to delete pQueryEncoder because ownership of it has
  // passed to *pQuery, so it is deleted when *pQuery is deleted.
  delete pQueryEncoder;
  delete pHashTable;
//  delete pSubjectMode;
//  delete pQueryMode;
  delete pQuery;
  if (pTask!=pPrintTask) delete pTask;
  delete pPrintTask;

} // ~processQuery

SequenceReader* findSequenceReader
( const string& fileName, const string& fileType, SequenceEncoder* encoder )
{

   if ( fileType == "fasta" )
   {
     return (SequenceReader*) new SequenceReaderFasta
       (fileName.c_str(),encoder,cerr);
   }
   else if ( fileType == "fastq" )
   {
     return (SequenceReader*) new SequenceReaderFastq
       (fileName.c_str(),encoder,cerr);
   }
   else if ( fileType == "" )
   {
     SequenceReaderMulti* pMulti = new SequenceReaderMulti(cerr);
     fillSequenceReaderMulti(*pMulti,encoder,fileName);
     if ( pMulti->getNumReaders() == 0 )
     {
       cerr << "Error: no recognized sources of sequence data were found in "
            << fileName << ".\n";
       exit (1);
     } // ~if
     return (SequenceReader*) pMulti;
   }
   else
   {
     cerr << "Error: " << fileType << " is not a recognized file type.\n"; 
     exit(1);
   }
   
} // ~findSequenceReader


// Name:      fillSequenceReaderMulti
// Arguments: SequenceReaderMulti& (out), string& (in), string& (in
// This uses the filename suffix to deduce what file type to create
// if fileName = *.fasta it creates an instance of SequenceReaderFasta and
// adds it to m.
// (Other file formats are TBD)
// If filename suffix does not exist or is not recognized, it assumes
// the filename must be a directory. It then calls itself for each filename
// found in the directory. By calling itself recursively in this way, all
// files in a directory tree are incorporated into m.
void fillSequenceReaderMulti( SequenceReaderMulti& m,
                              SequenceEncoder* encoder,
                              const string& fileName, 
                              const string& dirName )
{

  string::size_type dotPos(fileName.find_last_of('.'));
  string suffix
    ( (dotPos==string::npos)?"":fileName.substr(dotPos));  
  SequenceReader* pSeq(NULL);

  string fullPathName( (dirName=="")?fileName:dirName+"/"+fileName);
  //  cout << fullPathName << endl;

  if ( suffix == ".fail" )
  {
    return; // fail files picked up when their associated seq file is read
  }
  else if ( ( suffix == ".fasta" ) || ( suffix == ".fa" ) ) 
  { 
    cerr << "Creating SequenceReaderFasta from " << fileName << endl; 
    pSeq = new SequenceReaderFasta(fullPathName.c_str(),encoder,cerr);
    //    m.addReader(pSeq);
    //    m.addReader(new SequenceReaderFasta(fullPathName.c_str(),encoder,cerr));
  } // ~if
  else if ( suffix == ".fastq" ) 
  { 
    cerr << "Creating SequenceReaderFastq from " << fileName << endl; 
    pSeq = new SequenceReaderFastq(fullPathName.c_str(),encoder,cerr);
    //m.addReader(new SequenceReaderFastq(fullPathName.c_str(),encoder,cerr));
  } // ~else if
  // ... else ifs for all other formats ...

  if (pSeq!=NULL)
  {
    string filterFileName
      = ( (dirName=="")?"":dirName+"/")
      + ( (dotPos==string::npos)?fileName:fileName.substr(0,dotPos))
      + ".fail"; 

    ifstream* pFilterFile = new ifstream(filterFileName.c_str());

    if (!pFilterFile->fail())
    {
      SequenceReaderFilter* pFilter
        = new SequenceReaderFilter( pSeq, pFilterFile, cerr );
      pSeq=pFilter;
    }
    else delete pFilterFile;
    m.addReader(pSeq);

  }
  else // presume fullPathName is a directory, try to expand it
  {
    cerr << "Trying to open directory " << fullPathName  << " ..." << endl;
    DIR* pDir;
    if ( ! ( pDir = opendir( fullPathName.c_str() ) ) )
    {
      cerr << " ... failed! Exiting function." << endl;
      return; 
    } // ~if
    dirent* dirEntry;
    string entryName;
    while( dirEntry = readdir(pDir) )
    {
      entryName = (string) dirEntry->d_name;
      if ((entryName == ".")||(entryName=="..")) continue;
      // function calls itself recursively ... 
      fillSequenceReaderMulti(m,encoder,entryName,fullPathName);
    } // ~while
    closedir( pDir );
  } // ~else

  return;

} // ~findSequenceReader 

// End of file SSAHAMain.cpp

---