Global/SSAHAMain.h

Go to the documentation of this file.

/*  Last edited: May 30 13:49 2002 (ac2) */

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

// SSAHA : Sequence Search and Alignment by Hashing Algorithm
// Version 3.2, released 1st March 2004
// Copyright (c) Genome Research 2002

// SSAHA is free software; you can redistribute it and/or modify 
// it under the terms of version 2 of the GNU General Public Licence
// as published by the Free Software Foundation.
 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public Licence for more details.
 
// You should have received a copy of the GNU General Public Licence
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
// or see the on-line version at http://www.gnu.org/copyleft/gpl.txt

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

// Module Name  : SSAHAMain
// File Name    : SSAHAMain.h
// Language     : C++
// Module Author: Anthony J. Cox (ac2@sanger.ac.uk)

// Include guard:
#ifndef INCLUDED_SSAHAMain
#define INCLUDED_SSAHAMain

using namespace std;
// Description:

// Includes:

// NB it is good practise for #include statements in header files to be
// replaced by forward declarations if at all possible
#include <string>
#include <iostream>
#include <memory>
class SequenceReaderMulti;
class SequenceReader;
class HashTable;

// ### Class Declarations ###

struct QueryParameterStruct
{
  string queryName;
  string queryType;
  string queryFormat;
  string subjectName;
  string subjectType;
  string subjectFormat;
  int    queryStart;
  int    queryEnd;
  int    wordLength;
  int    stepLength;
  int    maxStore; 
  int    minPrint;
  int    maxGap;
  int    maxInsert;
  int    numRepeats;
  int    sortMatches;
  int    doAlignment;
  int    substituteThreshold;
  int    bandExtension;
  string queryReplace;
  string subjectReplace;
  //  string reportMode;
  string logMode;
  string saveName;
  bool   parserFriendly;
  bool   packHits;
  bool   reverseQuery;
  bool   runQuery;
  bool   printHashStats;
  friend ostream& operator<<(ostream& os, const QueryParameterStruct& qps )
  {
    return os << "queryName:\t" << qps.queryName   
       << "\nqueryType:\t" << qps.queryType   
       << "\nqueryFormat:\t" << qps.queryFormat   
       << "\nsubjectName:\t" << qps.subjectName 
       << "\nsubjectType:\t" << qps.subjectType  
       << "\nsubjectFormat:\t" << qps.subjectFormat  
       << "\nqueryStart:\t" << qps.queryStart  
       << "\nqueryEnd:\t" << qps.queryEnd  
       << "\nwordLength:\t" << qps.wordLength  
       << "\nstepLength:\t" << qps.stepLength   
       << "\nmaxStore:\t" << qps.maxStore    
       << "\nminPrint:\t" << qps.minPrint     
       << "\nmaxGap:  \t" << qps.maxGap     
       << "\nmaxInsert:\t" << qps.maxInsert     
       << "\nnumRepeats:\t" << qps.numRepeats
       << "\nsortMatches:\t" << qps.sortMatches
       << "\ndoAlignment:\t" << qps.doAlignment
       << "\nsubstituteThreshold:\t" << qps.substituteThreshold
       << "\nbandExtension:\t" << qps.bandExtension
      //       << "\nreportMode:\t" << qps.reportMode  
       << "\nqueryReplace:\t" << qps.queryReplace  
       << "\nsubjectReplace:\t" << qps.subjectReplace  
       << "\nlogMode:\t" << qps.logMode  
       << "\nsaveName:\t" << qps.saveName    
       << "\nparserFriendly:\t" 
       << (( qps.parserFriendly ) ? (string)"true" : (string)"false") 
       << "\npackHits:\t" 
       << (( qps.packHits ) ? (string)"true" : (string)"false") 
       << "\nreverseQuery:\t" 
       << (( qps.reverseQuery ) ? (string)"true" : (string)"false") 
       << "\nrunQuery:\t" 
       << (( qps.runQuery ) ? (string)"true" : (string)"false") 
       << "\nprintHashStats:\t" 
       << (( qps.printHashStats ) ? (string)"true" : (string)"false") 
       << "\n";
  } // ~operator<<

};

QueryParameterStruct defaultParams = 
{
  "",        // string queryName;
  "DNA",     // string queryType;
  "",        // string queryFormat;
  "",        // string subjectName;
  "DNA",     // string subjectType;
  "",        // string subjectFormat;
  1,         // int    queryStart;
  -1,        // int    queryEnd; 
  -1,        // int    wordLength;
  -1,        // int    stepLength;
  100000,    // int    maxToStore; 
  1,         // int    minToPrint;
  -1,        // int    maxGap;
  0,         // int    maxInsert;
  0,         // int    numRepeats;
  0,         // int    sortMatches;
  60,         // int    doAlignment;
  0,         // int    substituteThreshold;
  0,         // int    bandExtension;
  //  "tag",  // string reportMode;
  "default",  // string queryReplace;
  "tag",  // string subjectReplace;
  "cerr",    // string logMode;
  "",        // string saveName;
  false,     // bool   parserFriendly;
  false,     // bool   packHits;
  true,     // bool  reverseQuery; // made true the default TC 11.12.2001 
  true,      // bool   runQuery;
  false      // bool   printHashStats;
};  

// Class Name : CommandLineArg
// Description: This class and its subclasses interpret command line arguments.
// The interpretation is done by the static function parseCommandLine.
// It creates an instance of CommandLineArg for each valid command line option.
// Each instance recognizes a long and short option name and links to
// an int, string or bool. 
class CommandLineArg
{
  public:

  // 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 parameter 
  // values are sensible.
  static void parseCommandLine
  ( int numArgs, char* args[], QueryParameterStruct& queryParams );


  CommandLineArg( const string& nameLong, const string& nameShort ) :
    nameLong_( nameLong ), nameShort_( nameShort ) {}
  // Is the current argument equal to 'my' argument name?
  virtual bool isThisMe( const string& argName )
  { 
    return ( ( argName == nameLong_ ) || ( argName == nameShort_ ) );
  }
  virtual void addValue( const string& value ) = 0;
  protected:
  string nameLong_;
  string nameShort_;
};

// Class Name : CommandLineArgString
// Description: Recognize the name of an option that takes a string value, 
// and copy that value to the destination string
class CommandLineArgString : public CommandLineArg
{
public:
  CommandLineArgString
  ( const string& nameLong, const string& nameShort, string& destination ) :
    CommandLineArg( nameLong, nameShort ), destination_( destination ) {}
  void addValue( const string& value )
  {
    destination_ = value;
  }
protected:
  string& destination_;
};

// Class Name : CommandLineArgInt
// Description: Recognize the name of an option that takes an integer value, 
// and copy that value to the destination int
class CommandLineArgInt : public CommandLineArg
{
public:
  CommandLineArgInt
  ( const string& nameLong, const string& nameShort, int& destination ) :
    CommandLineArg( nameLong, nameShort ), destination_( destination ) {}
  void addValue( const string& value )
  {
    destination_ = atoi(value.c_str());
  }
protected:
  int& destination_;
};

// Class Name : CommandLineArgBool
// Description: Recognize the name of a boolean option. If recognized, set
// the destination bool to 'true'
class CommandLineArgBool : public CommandLineArg
{
public:
  CommandLineArgBool
  ( const string& nameLong, const string& nameShort, bool& destination ) :
    CommandLineArg( nameLong, nameShort ), destination_( destination ) {}

  virtual bool isThisMe( const string& argName )
  { 
    if ( ( argName == nameLong_ ) || ( argName == nameShort_ ) )
    {
      destination_ = true;
      return true;
    } // ~if
    return false;
  } // ~grabNextArg

  void addValue( const string& value )
  {
    destination_ = atoi(value.c_str());
  }
protected:
  bool& destination_;
};

// ### Function Declarations ###

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

class SequenceEncoder;

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



// 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 = "" );

// 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 );

static const char* helpText[] ={
"NAME\n",
" ssaha - performs rapid searching of DNA and protein databases\n",
"SYNOPSIS\n",
" ssaha [-h] [-help]\n", 
" - print this help message.\n",
" ssaha queryFile subjectFile [ -optionName_1 [optionValue_1] ] ...",
"   [ -optionName_n [optionValue_n] ] \n",
" - create a hash table from the sequences in subjectFile and", 
"   use it to search subjectFile for the sequences in queryFile.\n", 
" ssaha subjectFile [ -optionName_1 [optionValue_1] ] ...",
"   [ -optionName_n [optionValue_n] ] ...\n",
" - just create a hash table from the sequences in subjectFile.",
"   The -saveName option must be set (see OPTIONS).\n",
"DESCRIPTION\n",
" ssaha is a tool for rapidly finding near exact matches in DNA or protein",
" databases. The name is an acronym standing for Sequence Search and Alignment",
" by Hashing Algorithm. It works by converting a sequence database into a",
" hash table. This is then rapidly quizzed for hits, which are concatenated",
" into matches.\n",
"OPTIONS\n",
" Options may be specified by either their full names or short names and may",
" appear on the command line in any order.\n",
" Full Name      Short      Description\n",
"-queryFormat   -qf        Acceptable values:",
"                          fasta - fasta file",
"                          fastq - fastq file\n",
"                          Default value:",                        
"                          If not specified, attempts to deduce file type",
"                          based on the filename suffix as follows:\n",
"                          File suffix    Deduced file type",
"                          .fasta, .fa    fasta file",
"                          .fastq         fastq file",
"                          else assumed to be a directory of files, each of",
"                          whose names indicates the file type as specified",
"                          by the above rules.\n",
"-subjectFormat -sf        Acceptable values:",
"                          fasta - fasta file",
"                          fastq - fastq file",
"                          hash  - precomputed hash table\n",
"                          Default value:",                       
"                          If not specified, attempts to deduce file type",
"                          based on the filename suffix as follows:\n",
"                          File suffix    Deduced file type",
"                          .fasta, .fa    fasta file",
"                          .fastq         fastq file",
"                          else assumed to be a directory of files, each of",
"                          whose names indicates the file type as specified",
"                          by the above rules.\n",
"                          Note:",
"                          If -sf is set to hash, the -wl, -sl, and -ph",
"                          options will, if present, be ignored.\n", 
"-queryType     -qt        Acceptable values:",
"                          DNA",
"                          protein\n",
"                          Default value:",
"                          DNA\n",
"-subjectType   -st        Acceptable values:",
"                          DNA (in which case queryType must also be DNA)",
"                          protein",
"                          codon (i.e. do 6 way DNA to protein translation)\n",
"                          Default value:",
"                          DNA\n",
"-hashStats     -hs        Show information about the hash table currently",
"                          in use.\n",
"-parserFriendly -pf       Show one match per line as a set of tab delimited",
" (a.k.a. perlFriendly)    fields:\n",
"                          match direction: F forward, R reverse",
"                          query name",
"                          query start",
"                          query end",
"                          subject name",
"                          subject start",
"                          subject end",
"                          number of matching bases",
"                          percentage identity\n",
"-logMode       -lm        Controls the output of log information",
"                          Acceptable values:",
"                          cerr - send to standard error",
"                          cout - send to standard output",
"                          null - suppress log output",
"                          any other value sends log information to a file",
"                          of the same name\n",
"                          Default value:",
"                          cerr\n",
"-packHits      -ph        Store position of each word in a \"packed\"", 
"                          format comprising 32 bits per word. This halves",
"                          the size of the .body file at the expense of a",
"                          slight decrease in search speed.\n",
"-wordLength    -wl        Size in base pairs of the words used to form",
"                          the hash table. May vary from 1 to (assuming",
"                          sufficient RAM is available) 16.",
"                          Default value is 10.\n",
"-maxGap        -mg        Maximum gap allowed between successive hits for",
"                          them to count as part of the same match.",
"                          Default value is 0.\n",
"-maxInsert     -mi        Maximum number of insertions/deletions allowed",
"                          between successive hits for them to count as part",
"                          of the same match.",
"                          Default value is 0.\n",
"-maxStore      -ms        Largest number of times that a word may occur in",
"                          the hash table for it to be used for matching",
"                          expressed as a multiple of the number of",
"                          occurrences per word that would be expected",
"                          for a random database of the same size as the",
"                          subject database.",
"                          Default value is 10000.\n",
"-numRepeats    -nr        Maximum size of tandem repeating motif that can be",
"                          detected in the query sequence. This option may",
"                          produce faster and better matches when dealing",
"                          with data containing tandem repeats.",
"                          Defaults to 0, and must be less than or equal to",
"                          the word length.",
"                          Notes:",
"                          1. This option does nothing if -ph is also set.",
"                          2. To get the best results with this option, set",
"                          -mg to be at least equal to the word length.",
"                          Setting the -mi option may also help.\n",
"-minPrint      -mp        The minimum number of matching bases or residues",
"                          that must be found in the query and subject",
"                          sequences before they are considered as a match",
"                          and thus printed.",
"                          Default value is 1.\n",
"-queryStart    -qs        Specifies the number of the first query sequence to",
"                          be matched with the subject sequences (numbering of",
"                          both the query and subject sequences starts at 1).",
"                          Default value is 1.\n",
"-queryEnd      -qe        Specifies the number of the last query sequence to",
"                          be matched with the subject sequences. If not",
"                          specified, continues until the end of the query",
"                          sequence data is reached.\n",
"-reportMode    -rm        Specifies behaviour upon encountering unexpected",
"                          alphanumeric characters in query or subject ",
"                          sequences:\n",
"                          ignore - do nothing",
"                          report - report to standard error",
"                          replaceA   - silently replace character with 'A'",
"                          replaceG, replaceC, replaceT - as for replaceA",
"                          rrepA   - replace character with 'A' and report",
"                          rrepG, rrepC, rrepT - as for rrepA",
"                          Default value is `ignore.'\n",
"                          NB FOR VERSION 3.0, THE -reportMode OPTION HAS BEEN",
"                          SUPERCEDED BY THE -queryReplace AND -subjectReplace", 
"                          OPTIONS - SEE THE APPROPRIATE HELP ENTRIES\n",
"-reverseQuery  -rq        When matching the reverse strand of a query,",
"                          convert the positions of any matches found",
"                          into the coordinate frame of the forward strand.",
"                          Has no effect if queryType is set to protein.\n",
"-saveName      -sn        Specifies that the hash table must be saved before",
"                          the program exits. This option must be followed by",
"                          a string fileNameRoot. The hash table data is",
"                          saved into the files",
"                            fileNameRoot.head",
"                            fileNameRoot.body",
"                            fileNameRoot.name",
"                            fileNameRoot.size\n",
"                          Notes:",
"                          1. If no query file is specified (usage (iii)",
"                          above) it is an error not to set this option.",
"                          2. It is an error to set this option if", 
"                          subjectType is set to `hash.'",
"                          3. If the -ph option is also set, the -sn option",
"                          also produces a fileNameRoot.start file.\n",
"-sortMatches   -sm        Output only the top n matches for each query,",
"                          sorted by number of matching bases, then by",
"                          subject name, then by start position in the",
"                          query sequence.",
"                          Default value is zero, which outputs all matches",
"                          for each query and does no sorting.\n",
"-stepLength    -sl        Number of base pairs gap between words used to ",
"                          produce hash table. Ignored if a precomputed ",
"                          hash table is being used. Default value is ",
"                          equal to wordLength.\n",
"-queryReplace  -qr        Specifies behaviour upon encountering unexpected",
"                          alphanumeric characters in query sequences:\n",
"                          ignore - do nothing",
"                          report - report to standard error",
"                          A,G, etc. - replace with that character:", 
"                          must be A, G, C, T for DNA, or a valid IUPAC",
"                          amino acid code for protein.",
"                          Default: replace with 'A' for DNA, 'X' for protein\n",
"-subjectReplace -sr       Specifies behaviour upon encountering unexpected",
"                          alphanumeric characters in subject sequences:\n",
"                          ignore - do nothing",
"                          report - report to standard error",
"                          A,G, etc. - replace with that character", 
"                          Must be A, G, C, T for DNA or a valid IUPAC",
"                          amino acid code for protein",
"                          tag - `tag' the word so that it is not put", 
"                          into the hash table.",
"                          Default: tag\n",
"-substituteWords -sw      Look for single base/amino mismatches in words",
"                          that occur less than this many times more often",
"                          than would be expected for a random database of",
"                          the same size as the subject database.\n", 
"                          Only looks for:",
"                          purine (G-A)/pyrimidine (T-C) mismatches for DNA",
"                          mismatches with positive BLOSUM score for protein",
"                          Set to zero to switch this feature off.",
"                          Default value: 0 (switched off)\n",
"-doAlignment   -da        Produce a graphical alignment of the matching region",
"                          using banded dynamic programming. The alignment",
"                          will be formatted to the specified number of columns.",
"                          Set to zero to suppress alignments, otherwise",
"                          must be at least 20.",
"                          Default value: 80\n",
"-bandExtension -be        Specify size of the band to use for banded dynamic",
"                          programming, when producing a graphical alignment.",
"                          0 - diagonal only",
"                          n - n cells each side of diagonal",
"                          Only has an effect when -be is nonzero",
"                          Default value: 0 (diagonal only)",
"OUTPUT FORMAT\n",
"When full alignments are requested (-da set to nonzero) the software produces",
"a line of information then a graphical alignment for each match found:\n",
"i) DNA against DNA (untranslated)\n",
"RF      p1_1a788a06.q1c 515     538     p1_1a788f11.p1c 493     516     24",
"100.00",
"Alignment score: 13",
"Q:000000515 tttt-tgagacggagtctcgctct",
"            ||||x||||||xx|||||||||||",
"S:000000493 ttttttgagacaaagtctcgctct\n",
"ii) Protein against protein\n",
"FF      SW:PPSA_AERPE   625     642     SW:PPSA_METTH   577     592     8  ",
"44.44",
"Alignment score: 31",
"Q:000000625 KGGEKYETLDERNPMIGW",
"            x|||x |xx |x|||x||",
"S:000000577 EGGEN-EPY-EHNPMLGW\n",
"iii) Protein query against translated DNA subject ",
"(format for translated DNA query against protein subject is similar)\n",
"FR      SW:PPS2_HUMAN   467     473     p1_1a788c10.q1c 507     529     7",
"100.00",
"Alignment score: 22",
"Q:000000467 E..E..G..--V..L..D..P..",
"            |||||||||  ||||||Nxx|||",
"S:000000507 gaggagggcatgtattaaaccca\n",
"iv) DNA against DNA (translated)\n",
"FF      p1_1a788a01.p1c 52      76      p1_1a788b11.p1c 314     336     24",    
"96.00",
"Alignment score: 3",
"Q:000000052 atggtatgtctttcttttact-agat-",
"            ||xV..C..L..S..F..    R..  ",
"S:000000314 attgtttgtctctccttc---taga-a\n",
" From left to right the fields in the match information line are as follows:\n",
" i)    First character: query direction (F forward, R reverse)",
"      Second character: subject direction (F forward, R reverse)",
" ii)   query name",
" iii)  query start",
" iv)   query end",
" v)    subject name",
" vi)   subject start",
" vii)  subject end",
" viii) estimated number of matching bases",
" ix)   estimated percentage identity\n",
" The last two quantities are not exact values, they are approximations used",
" to order the matches (if requested) before the full alignment is done.\n",
" With the alignments switched off (-da 0), an entry like the one below is", 
" produced for each sequence in the query.\n",
" Matches For Query 6 (653 bases): p1_1a788a03.q1c\n", 
" F 6 : p1_1a788a03.q1c  Bases: 650   Q: 1 to 650    S: 1 to 650     100.00%",
" R 5 : p1_1a788a03.p1c  Bases: 100   Q: 22 to 121   S: 501 to 600   100.00%\n",
" The top line shows the query number, name and size of the query sequence.",
" Below that is one line for each match found in the subject database. From",
" left to right, the entries on these lines are as follows:\n",
" match direction: F forward, R reverse",
" subject number",
" subject name",
" number of matching bases",
" query start",
" query end",
" subject start",
" subject end",
" percentage identity\n",
"Notes:\n",
" 1. The output format is different if the program is run with the", 
" -parserFriendly option set. See the description of that option for details.\n",
" 2. Because SSAHA works by looking for whole-word matches, the `number of",
" matching bases' and `percentage identity' fields must be considered as lower",
" bounds on the true values of these quantities.\n", 
"FURTHER INFORMATION\n",
" The SSAHA home page is at http://www.sanger.ac.uk/Software/analysis/SSAHA/\n",
" Zemin Ning, Anthony. J. Cox and James C. Mullikin. SSAHA: A Fast Search",
" Method for Large DNA Databases. Submitted to Genome Research.\n",
"@"
};





// End of include guard:
#endif

// End of file SSAHAMain.h

---