/***************************************************************************** # Copyright (C) 1993-1996 by Phil Green. # All rights reserved. # # This software is part of a beta-test version of the swat/cross_match/phrap # package. It should not be redistributed or # used for any commercial purpose, including commercially funded # sequencing, without written permission from the author and the # University of Washington. # # This software is provided ``AS IS'' and any express or implied # warranties, including, but not limited to, the implied warranties of # merchantability and fitness for a particular purpose, are disclaimed. # In no event shall the authors or the University of Washington be # liable for any direct, indirect, incidental, special, exemplary, or # consequential damages (including, but not limited to, procurement of # substitute goods or services; loss of use, data, or profits; or # business interruption) however caused and on any theory of liability, # whether in contract, strict liability, or tort (including negligence # or otherwise) arising in any way out of the use of this software, even # if advised of the possibility of such damage. # #*****************************************************************************/ /* SWAT: a program for searching one or more DNA or protein query sequences, or a query profile, against a sequence database, using (an efficient implementation of) the Smith-Waterman or Needleman-Wunsch algorithms, with linear (affine) gap penalties. For each match an empirical measure of statistical significance derived from the observed score distribution is computed. 5/11/94 N.B. NW option will not work correctly until function call has been made similar to that for s-w. Usage: swat query_file_name database_file_name [-option value] [-option value] ... Available options are : option name value to be provided default value -matrix (or -M) name of file containing score matrix or profile BLOSUM50 -gap_init gap initiation penalty (= penalty for -12 first residue in a gap) -gap_ext gap extension penalty (= penalty for -2 each subsequent residue) -ins_gap_ext insertion gap extension penalty (= penalty for gap_ext each subsequent residue, for insertions in subject relative to query) -del_gap_ext deletion gap extension penalty (= penalty for gap_ext each subsequent residue, for deletions in subject relative to query) ) -end_gap gap extension penalty at ends of -1 sequences (N-W algorithm only) -E E-value cutoff; alignments with an E-value 1.0 higher than this will not be displayed. -z z-score cutoff; alignments with a z-value 6.0 lower than this will not be displayed. (See below for description of z-scores). -N maximum no. of displayed alignments 20 -truncatedb no. of entries at start of db that will be searched; useful for testing purposes only. The default is to search the entire database. -nw [flag; no value should be provided. Indicates that Needleman-Wunsch algorithm should be used instead of Smith-Waterman] -raw [flag indicating that z and E-values should not be computed, and matches should be ranked by their raw alignment scores.] -file [flag to create an output file for each query, giving the "complete" results of the search (but without alignments). This option is useful primarily for testing the relative power of different score matrices and gap penalties and will not be needed by most users. The file name is constructed automatically for each query and should NOT be given on the command line; it consists of the query sequence ID, followed by the suffix ".allscores". On completion of the run the file contains one line for each database entry giving the ID, length, score, and z-score. Entries are sorted by decreasing z. This file is produced independently of the usual program output and is unaffected by the -z, -N, or -E settings.] The query file (which may contain more than one query sequence, unless a profile search is to be performed in which case it can contain only one query sequence) and database file must each be in FASTA format: Each entry must have a single header line with the leading character '>' followed immediately by the ID of the sequence; the remainder of this line may contain an optional description. The header line is followed by lines containing the sequence. Any non-alphabetic characters are ignored, lower case letters are automatically converted to upper case, and any letter not in the score matrix alphabet is converted to the last character in this alphabet (typically an 'X' or '*'). N.B. NO LONGER CASE SENSITIVE -- all row and column labels should be capital letters. An appropriate matrix must be provided, in order to score alignments between the query (or profile) and each "subject" (database sequence). The matrix may be either a conventional score matrix, or a profile, but must contain only integer score values. The format of the matrix file is as follows: Any blank lines, and lines beginning with # are ignored. A line beginning with the word "GAP" is taken to contain the gap initiation and gap extension penalties (a pair of integers separated by spaces); this line must precede the line containing the column labels. The values given here will be overridden by command line parameters. The first non-blank line not beginning with # or GAP is assumed to contain the column labels, each of which should be a letter or '*'; there should be one column for each distinct type of residue that may appear in a subject sequence. The same label should not be used for two different columns. Succeeding lines should contain the rows of the matrix. Each row may optionally include a label (i.e. a single letter or '*' at the beginning of the line); if these are omitted, the row labels are taken to be the same as the column labels. If the matrix is a profile, row labels must be included, a single query sequence must be provided in the query file whose positions correspond 1-1 to the rows of the profile, and the sequence of this query must match the list of row labels (the only role of the query sequence in this case is to identify positions in printed alignments -- it plays no role in the analysis). If the matrix is a conventional score matrix, there should be one row for each type of residue that may appear in the query. Given an alignment, the score for an aligned pair of residues is determined by looking up the matrix entry whose column label matches the subject residue, and whose row corresponds to the query residue, either by position (if the matrix is a profile) or by label (if it is a conventional score matrix). Case is significant (so 'A' does not match 'a'). Residues for which there is no matching label are assigned to the last column (or row). Note that three different options, E, z, and N, influence how many alignments are displayed. If a single one of these is set on the command line it becomes the decisive criterion; if more than one is set, the most stringent is used. If none are set, the most stringent default value for the three options is used. By default, matches are sorted and displayed by decreasing z-score (see below), or equivalently E-value, rather than raw alignment score, unless the flags -raw or -nw are set. N.B. For Needleman-Wunsch searches, or if the database is small, z- and E-scores are not computed and the matches are instead ranked by the raw alignment score. Current restrictions: (i) In Smith-Waterman searches, if the alignment score for a database entry is less than -1 * gap_init it will be reported as 0. This restriction can be relaxed (at some cost in execution time) by moving the position of the corresponding block of statements in the source code. It is irrelevant for most purposes unless a very large gap_init value is used. (ii) ints must be at least 32 bits (so this program may not work with some compilers on DOS machines). (iii) The statistical analyses currently assume that "chance" scores are very unlikely to exceed 200. This assumption is reasonable when protein searches are done using the usual score matrices; it may not be reasonable for searches of nucleotide databases. It can easily be relaxed by changing appropriate program parameters. (iv) z- and E-values are only computed for Smith-Waterman searches. They may be unreliable for some combinations of score matrix / gap penalties. (see below), or if the database is small. (v) A single set of gap penalties is allowed; i.e. they cannot be made position-specific. Output: Mostly self-explanatory. includes a summary of the statistical analysis of the score distribution, and the alignments, scores and E-values of the best matches. Matches are ranked by decreasing z-score. Statistical procedures/ ranking of matches: (N.B. The following applies only to Smith-Waterman searches.) When the database sequences vary greatly in length (as is usually the case), theoretical arguments and empirical studies suggest that the raw alignment score is less sensitive, for discriminating true homologies from chance matches, than an adjusted score which has a null distribution independent of the length of the database sequence. SWAT attempts to compute such an adjusted score in the form of a "z-score", defined for a given database sequence by z = (S - f(n)) / sqrt(g(n)) where S and n are the raw alignment score and length of the database sequence, and f(n) and g(n) are the predicted mean score and variance for sequences of length n. z-scores thus should have mean 0 and standard deviation 1, independent of n. Like the raw score S, but unlike the E-value, z should be relatively independent of database size (although it may depend somewhat on the characteristics of database sequences); moreover it has a more easily interpretable scale than the raw score. At present, f and g are estimated from the observed score distribution for the database search, by regressing the observed scores and squared residuals, respectively, against log(n). This functional form has some theoretical justification (see below) and appears to work well for the score matrices and gap penalties in common use, although we are still exploring whether other functions may work better in some situations. (The regressions are actually performed twice, the results of the first one being used to identify high-scoring outliers (e.g. true homologies) which are then eliminated prior to performing the second regression). When the search has "local" behavior (i.e. the score distribution conforms to Karlin-Altschul type formulae), then known theoretical results imply that f(n) is in fact a linear function of log(n), g(n) is constant, and z is distributed as the (unique) Gumbel distribution with mean 0 and variance 1, independently of n. Consequently it is easy to convert z-values to E-values in this case, and they are monotonically related to each other. For searches that do not have local behavior, it is not known what the appropriate theoretical forms for f and g are, nor what the distribution of z is nor even whether this distribution is independent of n. Nonetheless z-scores should still be preferable to raw scores, and consequently the same procedures for computing z and E are used by SWAT in this case. They appear to work reasonably well at least for for PAM250 and gap penalties -12 and -2, a particular combination having non-local behavior that Pearson has found to be optimal for Smith-Waterman searches. Our tests using random queries indicate that the E-values obtained in this case tend to be somewhat conservative. A second procedure for estimating statistical significance is also included in the code but no longer used. It explictly assumes the search has local behavior and is based on fitting an extreme-value (Gumbel) distribution to the observed score distribution using maximum likelihood estimation (inspired by, but simpler than, Mott's procedure). (The idea of using the observed score distribution to estimate significance appears to be due to Collins and Coulson). The top part of the distribution (by default, the top 5% of the scores) is ignored, and a censored Gumbel distribution is fitted to the remainder to get estimates of lambda and K (in the notation of Karlin and Altschul). (This is different from Mott's procedure, which uses an outlier detection criterion that may miss many true positives. Also, we assume a simpler form for the distribution that does not allow lambda to vary for different database sequences; this makes it possible to use an iterative algorithm for finding the maximum likelihood parameter estimates that should converge much more rapidly.) For brief description of the implementation of the s-w algorithm, see header to smith_wat.c.